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ABSTRACT 



The Fermi-Pasta-Ulam (FPU) chains of particles in thermal equilibrium are stud- 
ied from both wave-interaction and particle-interaction points of view. It is shown 
that, even in a strongly nonlinear regime, the chain in thermal equilibrium can 
be effectively described by a system of weakly interacting renormalized nonlinear 
waves. These waves possess (i) the Rayleigh- Jeans distribution and (ii) zero cor- 
relations between waves, just as noninteracting free waves would. This renormal- 
ization is achieved through a set of canonical transformations. The renormalized 
linear dispersion of these renormalized waves is obtained and shown to be in excel- 
lent agreement with numerical experiments. Moreover, a dynamical interpretation 
of the renormalization of the dispersion relation is provided via a self-consistency, 
mean-field argument. It turns out that this renormalization arises mainly from the 
trivial resonant wave interactions, i.e., interactions with no momentum exchange. 
Furthermore, using a multiple time-scale, statistical averaging method, we show that 
the interactions of near-resonant waves give rise to the broadening of the resonance 
peaks in the frequency spectrum of renormalized modes. The theoretical prediction 
for the resonance width for the thermalized /9-FPU chain is found to be in very good 
agreement with its numerically measured value. Moreover, we show that the dy- 
namical scenario for thermahzed /3-FPU chains is spatially highly locahzed discrete 
breathers riding chaotically on spatially extended, renormalized waves. We present 
numerical evidence of existence of discrete breathers in thermal equilibrium. 
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0.1 Thesis outline 

The thesis is organized as follows. In Chapter [H we provide a historic reference 
to nonlinear science and the FPU problem in particular. In Chapter [21 we present a 
short overview of the theory of Hamiltonian mechanics. The methods of Hamiltonian 
mechanics will lead us to the dynamical description of the FPU chains. In Chap- 
ter [31 a brief introduction to the equilibrium statistical mechanics is given. There, 
we define the notion of microcanonical and canonical ensembles. In Chapter [H we 
present a formalism used in wave turbulence that we will later apply to the FPU 
system in order to give it a wave description. In Chapter [5l we describe symplectic 
integrators, which are numerical algorithms that conserve the symplectic structure 
of Hamiltonian systems. In Chapter [6l we give an introduction to chaos and dis- 
cuss a few examples that demonstrate the chaotic behavior of nonlinear systems. 
In Chapter [3, we rewrite the /3-FPU chain as an interacting four-wave Hamiltonian 
system. We demonstrate how to describe a strongly nonlinear system as a system 
of waves that resemble free waves in terms of the power spectrum and vanishing 
correlations between waves. We show how to construct the corresponding renor- 
malized variables with the renormalized linear dispersion. In Chapter [HI we study 
the dynamics of the chain numerically and find excellent agreement between the 
renormalized dispersion, obtained analytically (in Chapter [7]) and numerically. In 
Chapter [9l we describe the resonance manifold analytically and illustrate its control- 
ling role in long-time averaged dynamics using numerical simulation. In Chapter [TOl 
we derive an approximation for the renormalization factor for the linear dispersion 
using a self-consistency condition. In Chapter [Til we study the broadening effect of 
frequency peaks and predict analytically the form of the spatiotemporal spectrum 
for the /9-FPU chain. We also provide the comparison of our prediction with the 
numerical experiments. In Chapter [121 we discuss the energy localization in the 
/3-FPU chain in the form of discrete breathers. It was previously know that dis- 
crete breathers arise in the transient to thermal equilibrium under certain initial 
conditions. Here, we numerically demonstrate that discrete breathers also persist in 
thermal equilibrium. We present the conclusions in Chapter [131 
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CHAPTER 1 

Fermi-Pasta-Ulam problem as a part of nonlinear science 



In this chapter, we provide a general motivation for studying nonhnear science. In 
particular, we discuss the Fermi-Pasta-Ulam problem, its history, formulation, and 
various implication in physics and mathematics. 

1.1 Nonlinear macroscopic systems 

Study of nonlinear systems is a rich and fast growing direction in science [8]. 
Perhaps the key to the endless rich discoveries in nonlinear science arises from the 
fact that there are no general methods and approaches that can be used to uni- 
versally characterize any nonlinear system. Indeed, it is hard to imagine that a 
general algorithm can be developed to solve typical nonlinear equations. Since most 
of the phenomena that one finds in nature are nonlinear, one needs to look for 
special analytical methods of resolving particular nonlinear problems. For exam- 
ple, the simplest model of a physical system that everyone studies in school, the 
mathematical pendulum, is a nonlinear system. As we all know, in this case the 
linearization provides a very good approximation of the pendulum motion in the 
small-amplitude regime. Thus, one of the methods of treating a nonlinear system is 
to reduce it in some sense to a "close" linear system. However, the main drawback 
of this simplification is in losing all the rich nonlinear phenomena after linearizing. 
In addition, linearization does not always provide a physically meaningful approx- 
imation. Consider a water fiow in the pipe, when the velocity of the water is high 
and the flow is characterized as turbulent. Then we can hardly approximate the 
flow by a linear laminar water motion. Similarly, the formation of eddies and other 
spatially localized structures in the ocean can not be explained in the framework of 
spatially extended linear dispersive waves. 

Furthermore, many physical systems are not only nonlinear, they also have a 
large number of degrees of freedom. For example, the number of the air molecules in 
a room is of the order of Avagadro number, which is approximately equal to Na ~ 
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Figure 1.1: FPU chain. 

6.02 X 10^'^. Of the same order is the number of atoms in a crystal of a semiconductor. 
In biology, the DNA consists of millions of base pairs of molecules. In general, the 
study of many-body systems, or the so-called macroscopic systems, is a subject 
of statistical mechanics, which will be briefly introduced in the following chapter. 
All these and many other examples suggest that the study of many-body nonlinear 
systems can have important ramification in mathematics, physics, engineering, and 
biology. In order to characterize a behavior of a complicated physical system, one 
can first study its simplified model, which carries major features of its "parental" 
system while it is simple enough to be treated analytically or/and numerically. The 
celebrated Fermi-Pasta-Ulam (FPU) chain is one such system. 

1.2 History 

The study of discrete one-dimensional chains of particles with the nearest- 
neighbor interactions provides insight to the dynamics of various physical and bi- 
ological systems, such as crystals, wave systems, and biopolymers [10,48,57]. Let 
us first introduce the FPU model, which is going to be the main object of study 
in this thesis. Consider a one-dimensional chain of identical particles coupled with 
identical nonlinear springs as shown in Fig. II. 1[ This chain of oscillators was first 
introduced in the numerical experiment designed by physicist Enrico Fermi, com- 
puter scientist John Pasta, and mathematician Stan Ulam. Their last names gave 
the acronym FPU. The significance of the discovery made by FPU is manifested in 
a large number of theories and different branches of nonlinear science that appeared 
in the last fifty years. Here, we only present a concise description of the history 
of the problem, various attempts to resolve it, and some of its applications. More 
detailed information can be found in [20]. 
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In the early 1950s, the computer MANIAC I (Mathematical Analyzer Numer- 
ical Integrator And Calculator) was built and awaiting a significant question to re- 
solve. FPU proposed to use MANIAC I to integrate a one dimensional many-particle 
system, i.e., the FPU chain described above. The motivation for this experiment was 
to verify the fundamental beliefs of the statistical mechanics such as equipartition of 
energy among the degrees of freedom and ergodicity. Equipartition of energy means 
that all the degrees of freedom of the system contain the same amount of energy in 
average. And the ergodic hypothesis states that the time average of an observable 
is equivalent to the phase space average, i.e., average over phase space of the same 
nonlinear system taken at one time. The FPU chain with the third order nonlin- 
earity in the potential of the springs was simulated by a programmer named Mary 
Tsingou. The results were extremely surprising. Despite the expected equiparti- 
tion of energy, a very regular almost periodic behavior was observed. The system 
was initiated with only the first fundamental mode excited and after some time 
(~197 longest linear periods) the energy of the first fundamental mode recovered to 
within 3% of its initial value. This phenomenon was called "FPU recurrence" . In 
1955, the report of this experiment was distributed among the limited number of re- 
searches [18]. However, before this preprint was published, Fermi died. This sudden 
death prevented the preprint from being published — Pasta and Ulam could not 
publish the paper with Fermi's name on it since he had neither read or approved the 
manuscript. On the other hand, they could not publish it without his name since 
he was one of the creators of the whole idea. The report was published with the 
collected papers of Fermi only a decade later. 

The paradoxical behavior of the FPU chain initiated various attempts to re- 
solve the problem. Since the system exhibited a ncar-integrable behavior, the at- 
tempts to find a "solution" were made via approximating the system by completely 
intcgrable systems. (Note that the system is called completely integrable if it can 
be described by a Hamiltonian that is a function of the momenta only [44] .) Mar- 
tin Kruskal and Norman Zabusky noticed that the continuum approximation of the 
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discrete FPU problem is a famous Korteweg-deVries (KdV) equation 

ut + uu^ + u^xx = 0. (1.1) 

See Appendix |X] for the details on how to obtain the KdV equation using the small 
amplitude, long wavelength approximation of the FPU chain. The KdV equation is 
completely integrable and, moreover, possesses special solutions called solitons [59]. 

Solitons were first discovered in water canal by Scottish engineer John Scott 
Russell in 1834. He was watching horses pull a barge along the Union Canal in 
Edinburgh when the rope to the barge broke. The barge suddenly dipped into the 
water, which created a stable wave that set off up the canal with very little change 
in shape. Russell rode his horse along the canal for several kilometers watching the 
wave's progress. However, the significance of solitons in physics was only understood 
after the solitons were found to be the solutions of PDE's that describe physical 
phenomena, such as the KdV equation. 

Now we return to the connection of the KdV equation and the FPU chain. 
The numerical integration of the KdV equation with periodic boundary condition 
and the initial condition in the form of the one cycle of the cosine (the first Fourier 
mode as in FPU experiment) revealed an interesting dynamical behavior. The cosine 
first transformed into a number of spatially localized pulses, solitons, which after 
some time superposed back to the initial cosine. Although, this approach is only an 
approximation and does not provide a rigorous treatment of the FPU phenomenon, it 
gives an intuitive explanation of the observed phenomenon. However, the discovery 
of soliton that is a pure example of a coherent structure was significant by itself. 
There are numerous examples of coherent structures in nature [8]: from the giant 
(~ 10* meters) Red Spot in the atmosphere of Jupiter to the microstructures (~ 
10^^ meters) in crystals. After the solitons were found in nature, the whole class 
of completely integrable nonlinear partial differential equations became a rapidly 
growing research subject in mathematical physics. 

Further studies of the FPU model revealed that higher strength of nonlinearity 
induces irregular dynamics as opposed to the ordered, almost integrable behavior 
observed in the initial experiment. Thus, the study of another characteristic of many 
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nonlinear systems, chaos, was also influenced signiflcantly as a result of the FPU 
discovery. The existence of a certain threshold of the total system energy, called 
stochasticity threshold, which roughly separates the regimes with near-integrable 
and chaotic dynamics, was studied both numerically and analytically [11, 25]. If the 
total energy of the system is below the stochasticity threshold, then the recurrent 
behavior is observed as in the FPU experiment. However, when the total system 
energy increases and takes values above the stochasticity threshold, the dynamics 
of the chain becomes chaotic, and eventually the total energy becomes equally dis- 
tributed among all the degrees of freedom as was expected in the FPU experiment. 
A number of questions still remain open. The transition from the near-integrable 
to chaotic regime is not well understood. In particular, the behavior of the stochas- 
ticity threshold is not known in the thermodynamic limit (i.e., when the number 
of degrees of freedom goes to infinity). Moreover, the route to thermahzation (and 
hence energy equipartition) in the chaotic regime is not fully characterized. 

The near-integrable behavior of the FPU system with weak nonlinearity is 
closely intertwined with the celebrated Kolmogorov- Arnold-Moser theorem [44] . The 
theorem essentially says that for a small perturbation of the non-degenerate inte- 
grable system most of the invariant tori survive. That is to say, a small perturbation 
of the non-degenerate integrable system exhibits a near-integrable behavior. If, how- 
ever, the strength of the perturbation is increased, the invariant tori do break and a 
transition to chaos is observed. The theorem was first formulated by A. Kolmogorov 
in 1954 (without a proof) and then proved independently by V. Arnold in 1963 (for 
analytic Hamiltonian systems) and by J. Moser in 1962 (for twist maps — the area 
preserving maps with a phase twist that is radius-dependent). 



CHAPTER 2 
Hamiltonian Mechanics of discrete systems 



In this chapter, we discuss the classical mechanical approach of studying the behavior 
of a system of particles using the dynamical properties of the system. A full and 
complete study of the classical mechanics can be found in [33,37]. Suppose the 
system consists of particles, which are described by coordinates qj and A^ 
momenta Pj, where each component is a (i-dimensional vector and j is an integer in 
the range from 1 to A^. In this thesis, we will focus on a one dimensional system. If 
the system does not interact with anything — an isolated system as we will discuss 
in Chapter [3] — then it can be fully characterized by a Hamiltonian H{p, q), which 
is the total energy of the system. The dynamics of such a system is governed by the 
following canonical equations of motion 



It is often desirable to study the system in some new variables (P, Q) different from 
the initial ones, {p,q). The transformation {p,q) {P,Q) is called canonical if it 
preserves the form of the canonical equations of motion, i.e., Eq. (12.11) . 

For the weakly nonlinear systems, it is often convenient to make a transfor- 
mation from the physical space to the Fourier space since in the Fourier space the 
system can be viewed as weakly interacting waves. The discrete Fourier transfor- 
mation is defined via 



The Fourier transform of the real data, e.g., real qj, has the following symmetry 




(2.1) 




(2.2) 



6 



7 



property 

QN-k = Ql- (2.3) 
Using Eq. fl2.3l) . we show that transformation fl2.2p is canonical. 



• 1 . 2nikj 1 dH 2^ikj 1 dH dPi 



e 



N \/iV ^ dpj ViV dPi dpj 

j,«=i '« 

Similarly, we can show that 

Pk = -1^ (2.5) 
dQl ^ ' 

Yet another convenient transformation, which is often made, is a transformation to 
the so called normal modes. This is a linear transformation given by the following 
formula 

where Uk is an arbitrary positive function. Next, we show that transformation (12.61) 
is canonical if and only if 

i^k = l^N-k- (2.7) 

We have 

dok 1 fdPk dQk\ If dH dH 



dt \ dt dt J V dQl dP* 

1 / /dH dal dH daN-k\ ( dH dal dH da^-k 

v/2ZJ^ \\d^,dQl ^ da^-k dQl ) ~ '"^''Vd^^dP* ^ da^-k dP* 
1 f dH ^ fu^ _ ^ fu^\ ^ dH / j ujN-k _ ^ t^fc \ \ 



^/2uJ^ \dal\ V 2 V 2/ dajy-k^ \ 2 ^/2uJ^ 
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Then the equation of motion takes the canonical form 

dH , , 

if and only if property (12. 7p is satisfied. Here we have presented a few formal facts 
from Hamiltonian mechanics. We will use these transformations for studying the 
dynamical properties of the FPU chains. 



CHAPTER 3 
Equilibrium statistical physics 



Statistical physics is a branch of physics that studies systems with a large number 
of degrees of freedom (macroscopic systems). In principle, Newton's laws of motion, 
provide a way of describing the dynamics of all the particles in a given system. 
However, this is practically impossible for most physical macroscopic systems (for 
example, gas in a room) due to an enormous number of equations that have to 
be solved. Even if it can be solved, the macroscopic behavior of these trajectories 
requires a different conceptual framework. Statistical physics provides such tools — 
it uses the probabilistic approach, which is applicable in the case of a large number 
of particles in a system. 

3.1 Isolated systems and subsystems 

If a given system docs not interact with any other systems, then it is referred 
to as an isolated system. From our everyday experience we know that an isolated 
system reaches equilibrium state, i.e., the state when there are no temporal changes 
in any characteristics of the system. For example, if we pour some boiling water 
in the thermos then add a few ice cubes there and close the thermos, after some 
time the isolated system water+ice will equilibrate: the ice will melt and all the 
water in the thermos will have a uniform temperature, i.e., thermal equilibrimn will 
be reached. And if instead of the boiling water it was hot tea, the tea substance 
will also be uniformly distributed if we wait long enough. This shows that not 
only thermal but also chemical equilibrium is achieved. As another example of 
equilibration, we consider air in a room. Suppose, some perfume is sprayed in one 
corner of the room and then the room is left closed for some time. Eventually, the 
air will mix and the perfume will be uniformly distributed over the whole room. In 
these examples the use of the notion of an isolated system is consistent with the 
following facts. In the first case, the thermos docs not allow any thermal or chemical 
exchanges between the system water+ice with the outer world. In the second case. 
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walls, doors, windows, ceiling, and floor isolate the air in the room from the air 
outside the room so that the essential part of the perfume stays inside the room for 
a long time. 

Another type of system that is studied in statistical mechanics, is a system 
that is considered a part of a much large system. In this case this smaller system is 
referred to as a subsystem and the larger system is referred to as a thermal bath. It 
is important to point out that a subsystem should have a large number of degrees of 
freedom by itself. This makes the resolution of the dynamical equations of motion 
for each individual particle of the subsystem practically impossible. Therefore, a 
statistical approach should be applied to describe a subsystem. As an example, we 
can again consider the above described system water+ice in a thermos. However, 
now suppose that the thermos is not ideal and allows slow heat exchange between 
its content and the outer world. Suppose that this thermos with water+ice is left 
in a large room with a freezing temperature. Here, water+ice is a subsystem and 
the cold air in the large room is a thermal bath and the thermos plays a role of the 
interface between the subsystem and its thermal bath. First, the system water+ice 
will thcrmalizc but then in a much longer time (which depends on the goodness 
of the thermos) the water inside the thermos will freeze. Eventually, the content 
of the thermos will reach the temperature of the air in the cold room, i.e., of the 
thermal bath. This example demonstrates that depending on the physical situation 
and the observation time scales, the system water+ice can be either regarded as 
an isolated system or as a subsystem of a much larger system. If the thermos is 
good at isolating its content from the outer world, then the system water+ice can 
be treated as isolated. However, if the thermos is not perfect and the observation 
time is much longer than the time necessary to exchange heat between the content 
of the thermos and the cold air in the room, then the system water+ice is treated 
as a subsystem, which interacts with the thermal bath. 

3.2 Liouville's theorem 

In order to provide a statistical description of the system in thermal equi- 
librium we introduce a notion of a phase space. Suppose a system has degrees 
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of freedom, then its state at any time is described by coordinates and N mo- 
menta, qj and pj, respectively, with an integer j = {1, A^}, and the phase space 
is 2iV-dimensional space of points {q,p)- Time evolution of the system produces 
a trajectory {q{t),p{t)) in the phase space. Since the motion of a large system is 
very complex due to various interactions among the particles, the trajectory of the 
system after a long enough time should "cover" all the possible states of the system 
in the phase space. Mathematically it can be formalized by introducing the proba- 
bility density function (pdf) of all the possible states in the phase space. Consider 
a small volume element in the phase that contains all the points {q,p) such that 
Qj ^ i^p^j + ^^j) Pj ^ {P^P^j + dpj) for some point {q^^p^). Denote dw as 
a probability for the system to be in this small volume in the phase space. The 
probability dw can be expressed in terms of its pdf p via 



Let us derive one of the fundamental properties of the pdf of an isolated system, 
Liouville's theorem, which states that the the pdf p is constant along the phase 
trajectories. Formally, this statement comes from the conservation of the "number" 
of points in the phase space 



dw — p{q ,p )dqdp. 



Here p is normalized, so that 




(3.1) 



9p 
dt 



+ div(pt') = 0, 



(3.2) 



where v is the "velocity". In our case, the velocity is given by v = {q,p)'^- In the 
stationary state, we have 
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Then Eq. f l3.2p becomes 



N r 



E 



|-(p^.) + |-(pp.: 



And finally, taking into account canonical equations of motion, Eq. (12. ip . we obtain 
that the probability density function stays constant along the trajectories 



dp 



0, 



(3.3) 



where the so called Lagrangian derivative is defined via 



d 
di 



d_ 
di 



N 



•A -A 



(3.4) 



A direct consequence of Liouville's theorem is that p{q,p) is an integral of motion, 
and its logarithm is an additive integral of motion. If the system consists of two 
subsystems with pdf 's pi and p2 then the pdf of the whole system pu is a product of 
individual pdf 's of the subsystems, if the subsystems can be regarded as statistically 
independent. Therefore, we have 



lnpi2 = Inpi + lnp2- 



(3.5) 



There are in general seven independent additive integrals of motion as we know from 
mechanics. They are energy, three components of momentum and three components 
of angular momentum. Suppose we observe our system from a coordinate axis 
that is attached to the center of mass of the system, i.e., in this coordinate the 
total momentum vanishes. By a similar trick, we can eliminate the total angular 
momentum. Then, the only integral of motion of an isolated system would be its 
energy E{q,p). In the phase space, the motion of the isolated system becomes 
restricted to the surface that is given by the equation 



E{q,p) = Eo, 
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where Eq is the value of the total system energy. In the thermal equilibrium, we 
assume that all the points of this energy surface are equally probable. Then the pdf 
takes form 

p = const ■ 6{E{q,p) — Eq). (3.6) 

The ensemble of identical systems that are described by the probability measure (13. 6p 
is called the microcanonical ensemble. 

3.3 Entropy and the second law of thermodynamics 

Consider a small subsystem of a large isolated system in equilibrium. It turns 
out that the pdf of the subsystem in thermal equilibrium can be obtained analyt- 
ically. In order to do this we have to introduce entropy and temperature of the 
system. These quantities are related to each other and we start with entropy. As 
was noted above, the only additive integral of motion of an isolated system is its 
total energy. Therefore, for the subsystems of this isolated system we have 

lnps = as + bE,{q,p), (3.7) 

from where we conclude that p is a function of energy only. Then the total energy of 
the subsystem will be concentrated around a constant value E and the fluctuations 
around E of the total energy of the subsystem will be very small. Denote AqAp 
to be the volume of the phase space that consists of the phase points with energies 
close to E 

p{E)AqAp = 1. (3.8) 

This volume AqAp corresponds to the "number" of points of the phase space that 
have total energy close to E and the trajectory of the system should most of the 
time be inside the volume AqAp. Entropy is defined via 



(3.9) 
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where h = 6.58211915(56) x 10^^^ eV-s is Planck's constant and the constant {2TTh)^ 
ensures that entropy is a dimensionless quantity. Actually this factor has a deeper 
physical meaning — it is the smallest phase volume of a pair of degrees of freedom. 
Intuitively, entropy is a measure of disorder in the system — the higher values it 
takes, the more phase points have the energies close to E. In other words, with the 
same total energy, the state of the system with higher entropy is more disordered, 
i.e., has larger phase space volume to be in, than the state with the lower entropy. 
Let us obtain the expression for the entropy via the pdf p. Substituting AqAp from 
Eq. (13.81) into the definition of entropy (13. 9p . we have 

S = - \n{{27rhfp{E)]. (3.10) 

Next, using Eq. (13. 7p we obtain 

S = - ln(27r;i)^ -{a + bE) = - ln{2nh)^ - (In p{E)) = -{\n[{27Th)^ p{E)]) , 

where (...) denotes ensemble averaging over the probability measure given by p 
Thus we have derived an equivalent definition of entropy 

5* = -{ln[{27rhfp]) = - J p \n[{2Tr h) ^ p]dqdp, (3.11) 

From all observations it is known that if an isolated system is not in equilibrium 
then it will approach equilibrium as time progresses. This general law of nature is 
formalized by the entropy law or the second law of thermodynamics: 
If an isolated system is not in equilibrium at some moment of time then at 
subsequent moments of time its entropy will most probably monotonically 
increase. We have to point out that the statement of the entropy law as it is given 
here is not in contradiction with the time reversibility of the dynamical equations 
of motion, since here we only talk about the most probable state of the system. 
However, the issue with time reversibility of the dynamical equations and time 
irreversibility of the experimentally observed property of entropy, i.e., that entropy 
never decreases (if we disregard micro-fluctuations) is a deeper open issue and we 
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will not pursue this question here. For our purpose, we will only need the fact that 
in thermal equilibrium entropy has its maximum. 

3.4 Temperature and Gibbs measure 

Let us introduce one of the most important thermodynamic quantities, i.e., 
temperature. Suppose an isolated system consists of two subsystems, which are 
in equilibrium with each other. Denote Ei and E2 to be total energies of each 
subsystem and 5*1 and S2 to be entropies of each subsystem. For the total energy 
and entropy we have 

E — El + E2, 

S = Si{Ei) + S2{E2). 

Since £■ is a constant and E2 — E — Ei, total entropy 5" is a function of Ei only. 
The necessary condition for S to have maximum in thermal equilibrium is 

dS dSi dS2 dE2 _ dSi dS2 _ {o 



dEi dEi dE2dEi dEi dE2 

This property of the entropy of the subsystems is easily generalizable to any number 
of subsystems. Therefore, we can define the temperature T as 

dE T ^ ^ 

The temperatures of the two systems in thermal equilibrium are equal 

Ti = T2. 



Now we obtain the pdf of the subsystem, which is equilibrium with its thermal 
bath. Suppose that the subsystem together with the thermal bath is an isolated 
system. Then the microcanonical distribution for this combined system is described 
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by 

dw = const ■ 5{E + E' - E^^'^) dqdpdq'dp', (3.14) 

where E, E', and E^o) are the energies of the subsystem, of the thermal bath, and 
the combined system, respectively. Then the probability density function for the 
subsystem becomes 

p = const J 5{E + E' - E(°)) dq'dp'. (3.15) 

The integrand in Eq. (13.151) depends only on E', therefore, we can change the inte- 
gration variables from q',p' to E' via 

dq'dp' ^ ^^dE'. (3.16) 

Note that AE' is the length of the energy interval that corresponds to the phase 
space volume given by Aq'Ap'. Using the entropy definition given in Eq. (13. 9p we 
obtain 



P 



/S' 
^5{E + E' - E^°^)dE'. (3.17) 



In equilibrium, the fluctuations of E' are very small and the pdf of the energy of the 
thermal bath as a function of E' has a sharp peak around its average value. And the 
width of this peak AE' is practically independent of the energy of the subsystem E. 
Therefore, we can use the value E' = E^^^ in AE' and after integration Eq. (13.171) 
becomes 

p = const ■ e^'\E/=E(o)-E- (3.18) 

Since E is small compared to E^^\ we use the Taylor expansion with small parameter 
E in 5'(EW - E) 

5'(^(o) _E)= 5'(eW) - E^^^^ffp-. (3.19) 
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After combining Eqs. 03.181) and (13.191) and using the definition of temperature (13.131) 
we obtain the following form for the pdf of the subsystem 

p(q,p) = Zexp(^-^^y (3.20) 

and the so called partition function Z is defined from the normalization condi- 
tion ^l\f 

exp(-^%^] dqdp. (3.21) 



Distribution given by Eq. (I3.20p was found by Gibbs in 1901 and is referred to 
as the Gibbs distribution or canonical distribution or more commonly Boltzman 
distribution. 

The explicit form of the pdf given by Eq. (I3.20p is very convenient when the 
average characteristics of the various dynamical quantities have to be computed. As 
a simple example, we consider an atom with mass m and Hamiltonian 

H = ±^p^ (3.22) 
According to Eq. (I3.20p . the probability measure has the form 

dw = ^ exp (^-^^P^^ dp. (3.23) 
The average value of the kinetic energy of the atom can be easily computed, 

2m/ 2m/exp (-2^p2) 2' 

This results is generalized in the equipartition theorem, which states that for the 
system in thermal equilibrium described by a quadratic Hamiltonian the following 
relationship holds 

^m^) = SnT, (3.25) 
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where Xj is a degree of freedom, i.e., one of Qj or pj. 

However it is necessary to point out that practically both micro canonical and 
canonical distributions become identical when the number of degrees of freedom goes 
to infinity. The only difference between these distributions arises when one computes 
the fluctuations of the total energy around its average value. In the canonical 
distribution, these fluctuations are non-zero. On the contrary, in the microcanonical 
distribution, they are zero by definition. And as the number of degrees of freedom 
grows, the fluctuations of the total energy decrease as ~ l/y/N [51]. Prom the 
practical point of view, the calculations using the canonical distribution are much 
easier mathematically. We will use both notions of the microcanonical and canonical 
distribution when we study the statistical behavior of the FPU system from the wave 
point of view and we will see the equivalence of both approaches. 



CHAPTER 4 
Wave Turbulence 



Wave turbulence theory (WT) studies a statistical state of a system of nonlinear 
dispersive waves that weakly interact with each other, and their dynamics is de- 
scribed statistically. WT has been used for almost eighty years to provide a statis- 
tical description of various physical systems. Peierls initiated the methods of WT 
in [47], in which the kinetic equation for phonons in solids was obtained. Among 
other examples of WT are ocean, atmosphere, plasmas and Bose-Einstein conden- 
sates [5,24,50,62]. Perhaps the key discoveries in WT are made in application to 
oceanography [60-62] by V.E. Zakharov et al. There, it was argued that systems of 
dispersive waves develop a Kolmogorov type of turbulence in non-equilibrium state 
as opposed to the thermalized state as it was studied before. Kolmogorov- Zakharov 
non- equilibrium spectra that predict cascades of various excitations, play a central 
role in the modern development of WT. These spectra arise in the systems that are 
driven away from equilibrium by forcing and damping. The non-equilibrium situa- 
tions are the main focus of WT. However, in this thesis we only discuss the thermal 
equilibrium state of the FPU system. Nevertheless, we will use the ideas and meth- 
ods that are commonly used in WT. Therefore, we provide a brief description of 
WT here. 

In the general setting of WT, dispersive waves are governed by a Hamiltonian, 
such as 

H = J tu klttkl"^ dk + - yT^'^a^a^amfls^ms ^^^^'^'^'^"5) (4-1) 

where ak{t) describes the evolution of the kth wave mode in time, Uk is a linear 
dispersion, and T^'^ is an interaction tensor coefficient, which is considered to be 
small in the case of the weak coupling. The formal procedure for obtaining the 
Hamiltonian of type (14. ip is given in Chapter [2] and a more comprehensive discussion 
is provided in [62]. In Chapter [9], we will derive a discrete form of Hamiltonian (14.11) 
for the /3-FPU chain [Eq. (19. 3p ]. In Eq. (14. ip . we consider only the fourth order 
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interactions among the waves. The /?-FPU system that we will study is of the same 
type with quartic potential interactions. WT aims to derive the kinetic equation for 
the power spectrum defined by 



where (...) stands for averaging over an ensemble of initial data. This is usually 
achieved by combining the statistical description of the wave field with its dy- 
namic description using the Hamiltonian (14. ip . 

It is assumed that the wave field is near-Gaussian if the nonlinear interac- 
tions are weak. The near-Gaussian assumption leads to the following approximations 
of the fourth and sixth order correlators 



These approximations are crucial in making the closure in the hierarchy of equations 
for each order of small parameter. 

Now we outline the main steps in deriving the four-wave kinetic equation. The 
dynamical evolution of ak{t) is described by the following equation 




(4.2) 




nuni{5% + 5%), (4.3) 
rifcn,n^(5^(457 + 5^5^) + 5^(5^57 + (4.4) 



^ak = — = uju^k + / T^sf^o.mO's^ms dldmds. 



(4.5) 



'k 



We compute the time evolution of the wave action 




(4.6) 



where 



kl 



(4.7) 



ms 



Now, we use the approximation of the fourth order correlator and in the zeroth 
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order we obtain 

hk = 2nkQ [ Ttlni dl. (4. 



In Eq. f l4.8p . the integrand is real and therefore the RHS is equal to zero. Therefore, 
in order to find the higher order contribution, we should not split the fourth order 
correlator and should consider its time evolution as a whole. We obtain the following 
dynamic equation 

+ {uk + uji-ujm- uJs)^ J^s = 2fl^^{nmns{nk + rii) - nkni(nm + ris)). (4.9) 

Now we use the notion of the separation of linear and nonlinear time scales. We sup- 
pose that the integral contribution of the fast oscillations vanishes as time increases, 
which yields 

jki ^ ^T^si^mUsink + rii) - nkfiiin,^ + n,)) 

Uk + UJl - UJm - i^s + 10 

where the small term i5 was added. Physically, the term i5 represents friction, 
which is always present in any realistic system. However, this theoretic treatment 
has a profound effect on wave statistical dynamics and it is consistent with the fact 
that the resonances dominate the long time dynamics. Next, we use the following 
equality 



S = -^^(^)- (4-11) 

Finally, we obtain the four-wave kinetic equation 

nfc = 27r J \Tl^f{nmns{nk + rii) - ^^^^(n^ + n,))S'^J{uj^^) dldmds, (4.12) 

where = Uk + uji — Um — ujg- We note that the four- wave resonance conditions 
arise here: only those quartets of waves that are on or very close to the resonance 
manifold provide effectively mix of energy through the four wave interactions. And 
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the resonance manifold is described by the set of wave numbers that satisfy 

{k + I = m + s, 
(4.13) 
LOk + iOi = Um + UJg. 

We will study the interaction of waves in the thermalized state of the FPU 
system also using near-Gaussian assumption and notion of separation of linear and 
nonlinear time-scales. We will also study the resonances manifold of the FPU chain, 
which is described by the discrete version of the Eq. fl4.13p 



CHAPTER 5 
Symplectic integrators for Hamiltonian systems 



With the growing power of computers, it becomes more feasible to numerically simu- 
late many-body systems with a large number of degrees of freedom. However, usual 
numerical algorithms such as Euler scheme or Runge-Kutta methods do not conserve 
the total energy of a Hamiltonian system for long time simulations. Hamiltonian 
systems possess symplectic structures, which we define below. Therefore, efficient 
and precise algorithms that capture this crucial characteristic of Hamiltonian sys- 
tem, are needed. Here, we describe a class of such algorithms called symplectic 
integrators and provide all the necessary formulas and parameters for one particular 
method, i.e., the sixth order Yoshida method [58]. 

Suppose, we study a system of particles. At any given time the position of 
the system is described by the set of d ■ N coordinates and d ■ N momenta, which 
we will denote as a pair {q,p) assuming that both components are d-dimensional 
vectors. In Hamiltonian mechanics, the evolution of the system is given by the 
Hamiltonian function H{p,q) via the canonical equations of motion (12. ip . These 
equations define a time fiow of the phase space — to find the position (grjPr) of the 
system at time r one can integrate Eq. (12.11) up to t = r with the initial conditions 
(go,Po) given at t = 0. By definition, the time fiow of the phase space is symplectic 
if it preserves the differential form 



In order to construct the symplectic integrator, we provide a formal description of 



the Hamiltonian fiow given by Eqs. (12. ip . Define z = {p, q) and the Poisson bracket 
{-, ■} as 



00 = dp A dq. 



(5.1) 



{f,9} 




(5.2) 
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Then Eqs. (12. ip can be written in the form 

z = {z,H{z)}. (5.3) 

By introducing the notation = {•,//(■)} for the differential operator, Eq. (15.31) 
can be rewritten as 

i = Dhz, (5.4) 

and its formal solution t = r is given by 

z{t) = e^^"z{0). (5.5) 

Since the total energy is conserved, we can write 

H{z{T))=Hizm. (5.6) 

Suppose, the total energy of the system is a sum of the kinetic energy that is a 
function of p only and a potential energy that is a function of q only 

H{q,p) = T{p) + V{q). (5.7) 

Then the corresponding differential operators are denoted as 



jj ^ dTd_ 
dp dq'' 

n — —9Vd_ 
dq ap' 



(5. 



and thus the formal solution (15.51) becomes 

z{t) = e"(^^+^^)z(0). (5.9) 

Since the differential operators Dt and Dy are non- commutative, exponential of the 
sum of two operators in Eq. ( 15. 9p is not equal to the product of exponentials of the 
individual components. Instead, for any non- commutative differential operators A 
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and B, the approximate relationship holds 

e^(^+^) =e^V^ + o(r2). (5.10) 

Our goal is to build the n-th order symplectic scheme, and in order to achieve this 
we generalize Eq. (15.101) and construct the expansion of the form 

^r(A+B) ^ ^c^rA^d^rB^C2rA^d^rB ^...^ ^c,rA^d,rB _^ o(r"+^), (5.11) 

where ci . . . Cfc and di . . .dk are real numbers. Now, the solution for z{t) is approxi- 
mated by 

z'{t)= ^JJe^^-^V^^^j z(0), (5.12) 

Note that Eq. (I5.12p provides a symplectic mapping in the phase space, i.e., because 
it consists of a series of elementary symplectic mappings e'^-''^^ and e'^^'^^ . More- 
over, using Taylor expansions of the operators (15. 8p up to the first order in r, the 
corresponding elementary mappings are explicitly computable 



qU) = g(i-i) + rc,^(j9(^-^)), (5.13) 
dV 



pyjj — p 



0) ^ (5.14) 



dq 



for j = 1 to j = k. Naturally, the question arises whether the expansion of the 
form (15. lip exists for any order n. It turns out that for every even order n there 
exists at least one set of exact coefficients ci . . . Cfc and di . . .d^ so that Eq. ( 15.12^ 
provides an order n approximate solution for the dynamical equation (15. 3p . The 
general approach of finding the coefficients Cj and dj is as follows. We expand the 
LHS of Eq. (15. lip in powers of r up to the order n. Then we equate the coefficients 
of the corresponding powers on both sides of Eq. (15. lip and obtain a system of 
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equations for Cj and dj. The resulting coefficients for n = 6 are given by [58] 



Cl 


= 0.392256805238780 


di 


= 0.784513610477560, 


C2 


= 0.510043411918458 


d2 


= 0.235573213359357, 


C3 


= -0.471053385409757 


ds 


= -1.177679984178870, 


C4 


= 0.068753168252518 


di 


= 1.315186320683906, 


C5 


= 0.068753168252518 


da 


= -1.177679984178870, 


C6 


= -0.471053385409757 


de 


= 0.235573213359357, 


C7 


= 0.510043411918458 


d-i 


= 0.784513610477560, 


C8 


= 0.392256805238780 


ds 


= 0.0, 



To summarize, here we presented the numerical algorithm for solving the Hamilto- 
nian equations of motion (12. ip using the symplectic algorithm given by Eqs. (15.131) 
and (15.141) . Note that strictly speaking a symplectic algorithm may not preserve 
energy as well as an explicit Runge-Kutta method for a fixed time step r. How- 
ever, it is superior to an explicit Runge-Kutta method since the system energy is 
bounded when computed by a symplectic algorithm in contrast to the unbounded 
energy when computed by an explicit Runge-Kutta method. This property of a 
symplectic algorithm becomes important when a long time (statistical) behavior of 
a Hamiltonian system is studied. 



CHAPTER 6 
Chaos in dynamical systems 



Here we present a few general facts from the theory of chaos, which provides a 
language to describe random, turbulent, irregular behavior of the dynamical systems. 
A systematic and comprehensive discussion of chaos can be found, e.g., in [44]. 
Chaos arises in many dynamical systems. It was first noticed by H. Poincare when 
he studied the motion of the system of three celestial bodies and realized that 
very complex trajectories could arise. Thereafter, such scientists as G. Birkoff, 
M. Cartwright and J. Littlewood, S. Smale and A. Kolmogorov and many others 
investigated various aspects of chaotic dynamics. However, the importance of chaos 
was really appreciated for understanding real complex physical phenomena only 
when the computational power became widely available. 

6.1 Logistic map 

We start with an example of discrete maps in which both regular and chaotic 
behaviors can be observed. The logistic map is given by the following nonlinear 
transformation 



and A is a controlling (bifurcation) parameter. The map has two fixed points xi — 
and X2 = 1 — 1/4A. For the parameter range A > 1/4, the fixed point xi is unstable. 
For / G (1/4, 3/4), the fixed point X2 is stable and for all other / it is unstable. Here, 
the fact that the fixed point is stable means that a trajectory that is initiated at 
any point except for the set of point of the Lebesgue measure zero will eventually 
converge to the fixed point. Then this fixed point is the attractor of the map. 



(6.1) 



where 



f{x)^A\x{l-x). 



(6.2) 
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However, the attractor of the map can be more comphcated than just a point. The 
attractor, i.e., the set of points to which the set of all possible initial points converges 
(or becomes attracted), can be a very complicated object such as a fractal, i.e., whose 
Hausdorf dimension [21] is not an integer (then it is called a strange attractor). We 
will observe the route of the logistic map to chaos with the change of the parameter 



In Fig. 16.11 we show the attractor of the logistic map (16.21) as a function of /. 
We observe that when / crosses the value I = 3/4, the system encounters a pitch-fork 
or period doubling bifurcation [21] — the single line that represents the fixed point 
solution splits into two lines that represent a period two orbit. Then the attractor 
becomes two points instead of one in the case of a fixed point. 

After / is increased further, the bifurcation happens again but now to each of 
the two branches. After the nth bifurcation, the length of the period is 2". As n 
goes to infinity, the period becomes infinitely large and the system becomes chaotic. 
Here, this happens at /qo ~ 0.89286... For most / > /qo we have a case of a strange 
attractor. However, periodic orbits still exist even for / > /qo- But most of the orbits 
are chaotic for I > l^o- 

In general, a chaotic orbit is characterized by the following properties 

• it is not periodic, 

• it is not attracting to a periodic limiting set, 

• it has sensitive dependence on the initial conditions. 

A convenient diagnostic of the sensitive dependence on the initial conditions (or 
small orbit perturbations) is the Lyapunov exponent. Lyapunov exponent gives 
an average exponential rate of divergence of the two trajectories with close initial 
conditions. If we denote dxn = Xn+i — Xn then /i is a Lyapunov exponent if 




(6.3) 
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Figure 6.1: Attractor of the logistic map as a function of the parameter I. 



Thus we come to the strict definition of a Lyapunov exponent 



h = hm — IuIAatI, 



(6.4) 



where 



An = f'ix^)f'ixN^,)...fix,) 



(6.5) 



We note that if h > 0, then two initially close trajectories will diverge with the rate 



h, i.e., the case h > indicates chaos. In Fig. I6.2[ we demonstrate the dependency of 
the Lyapunov exponent of the logistic map on the controlling parameter /. Note that 
the periodic orbits correspond to the dips {h < 0) in the plot and the bifurcation 
points correspond to the zero values of the Lyapunov exponents. For the continuous 
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Figure 6.2: Lyapunov exponent of the logistic map as a function of the parameter 



trajectories we have the following generalization of the Lyapunov exponent 



h= \n\M'{x)\dij{x) 



(6.6) 



where |M'(x)| is a jacobian of the transformation Xn+i = M{xn) in the continuous 
limit dxn and fi{x) is an invariant measure [44]. 

To summarize, we have considered an example of a discrete map. We have 
observed that by changing the controlling parameter /, the orbits of the logistic map 
exhibit a change of behavior from regular with a periodic limiting cycle to chaotic 
with a strange attractor. Using this example, the important notions of the theory 
of chaos such as attractor and Lyapunov exponent are introduced. 
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6.2 Lorenz model 

We continue our discussion of chaotic systems with a continuous autonomous 
dynamical system — the Lorenz model. The model was originally obtained by 
truncating the Fourier expansion of the Navier-Stokes equations. This model also 
exhibits a change of behavior from regular to chaotic as the controlling parameter 
changes. The model is given by the following equations 

X = -a{x-y), (6.7) 
y = (^r- z)x- y, (6.8) 
z = xy — bz, (6.9) 

where x, y, and z are real functions of time and a, r, and b are real positive param- 
eters. We follow [21] and choose o" = 10 and 6 = 8/3 for our numerical experiments. 
This leaves r to be the only parameter that is varied. 

For < r < 1, the origin is the only attracting fixed point with all the orbits 
approaching the origin. For r > 1, two stable fixed points arise 

Ci,2 = (± VK^-l),±A/K^-l),r-l). (6.10) 

On the other hand, the origin loses its stability at r = 1. In Fig. 16.31 we observe an 
orbit that converges to the fixed point C2 for r = 5. If we increase r further then 
we find that for r > 27.74..., the fixed points Ci^2 become unstable and the behavior 
of most orbits become chaotic. The form of the nonlinearity prevents the orbits to 
diverge to infinity. In Fig. 16. 4[ we observe a typical behavior of an orbit for r = 30. 
Since the fixed points Ci^2 are now unstable, the orbit spirals outward around one 
of the fixed points. Then it reaches the vicinity of the other fixed point and spirals 
out again and so on. The behavior of most of the orbits in this parameter range is 
chaotic. In order to see that we can consider some discrete subset of the continuous 
trajectory and show that this discrete subset actually produces a chaotic mapping. 
In particular, we consider a trajectory z{t) and construct a sequence Zn of the local 
maxima of z{t). Then this sequence z„ Zn+i is approximately a one- dimensional 
map function, which is shown in Fig. 16.51 Note that the slope of this resulting one- 
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Figure 6.3: Orbit converges to a stable point C2 for r = 5. 

dimensional map is greater than 1 (compare with a black line in Fig. 16.51) . Therefore, 
\dzn+i/dzn\ > 1 and the Lyapunov exponent h is positive which indicates chaos. To 
complete the discussion we just add that the windows of periodic motion appear if 
the parameter r is increased. 

This example demonstrates that a system can exhibit both regular and chaotic 
behaviors and by changing a control parameter one can achieve both. 



6.3 Sensitive dependence on initial conditions 

One of the characteristic properties of chaotic dynamical systems is the sensi- 
tive dependence on the initial conditions. Suppose we have a dynamical system 

^ = F[x{t)], (6.11) 
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Figure 6.4: Chaotic behavior of the an orbit for r = 30. 

where x{t) is a vector solution. Then if the system is chaotic, any two trajectories 
that are initiated at two very close points will exponentially diverge from each other. 
This can be formulated more rigorously in the following way. Consider two points 
Xi(0) and X2(0) such that |xi(0) — X2(0)| is small. Then if the system (16.111) is 
chaotic, then the difference \xi{t) — X2{t)\ will grow with time exponentially, i.e., 
positive Lyapunov exponent, where Xj{t) is the solution of Eq. (16. lip with the initial 
value Xj{0). 

Therefore, even a small perturbation in the initial condition yields a completely 
different solution of a dynamical system in a chaotic regime. This property has a 
direct impact on the numerical computations. We always introduce some roundoff 
and truncation errors in both defining the initial conditions and then in computing 
the trajectory itself. Therefore, for chaotic systems, we can never obtain a numerical 
solution, which is close to the exact solution to the equation with a given initial data. 
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Figure 6.5: Zn+i as a function of z„ for the Lorenz model. The black line with slope 
1 is shown for comparison. 



If we use slightly perturbed initial data or different numerical integration schemes 
in the computations, for any reasonable longtime trajectory, i.e., T = 0{h^^), we 
will obtain a completely different numerical solution. 

Let us demonstrate the sensitive dependence on the numerical noise using the 
/3-FPU system. First, we show that the /5-FPU system is indeed chaotic, i.e., the 
Lyapunov exponent is always positive. In Fig. 16. 6[ we present the comparison of the 
successive approximations to the Lyapunov exponents of the /3-FPU system for four 
different values of the nonlinearity strength f3 = {0.25,1,4,16} and same system 
energy H = 100. The algorithm for computing the Lyapunov exponent of the /?- 
FPU system is given in Appendix [Bl Note that although the Lyapunov exponents 
are positive for any strength of nonlinearity, larger values of the Lyapunov exponent 
corresponds to larger values of nonlinearity. Therefore, we expect more chaotic 
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Figure 6.6: Positive Lyapunov exponents. 



behavior for strongly nonlinear systems. In Fig. 16.71 we compare the trajectories of 
the /3-FPU system with = 128 particles, simulated with the same parameters as 
in Fig. 16.6^ but with different time steps. The red curves correspond to dt = 0.01 and 
the blue curves correspond to dt = 0.02. Different time steps imply that different 
numerical noises were introduced during computation, which should result in the 
divergence of these two trajectories due to the chaotic nature of the /3-FPU. We 
make at least two observations after examining Fig. 16. 7[ (i) The /3-FPU system is 
indeed chaotic since the two trajectories diverge resulting in the positive Lyapunov 
exponent, h > 0. (ii) For higher values of nonlinearity strength (3 the system is more 
chaotic, meaning that this divergence occurs faster since the Lyapunov exponent is 
larger for stronger nonlinearity. 




Figure 6.7: Time evolution of q2 for the /3-FPU system with the different value of 
the nonlinearity parameter j3. 
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6.4 Shadowing theorem 

Sensitive dependence on the initial condition and numerical noise seems to 
destroy any hope for solving the chaotic dynamical systems on the computer since 
any tiny noise (such as numerical errors) will eventually bring in a huge error into the 
solution very quickly. However, the cure to this problem comes from a rigorously 
proved shadowing theorem [45]. Suppose we have an exact solution x{t) and a 
numerical solution x{t) and as we have seen their difference grows exponentially with 
time due to roundoff errors although a;(0) = x{0). The shadowing theorem says that 
there exists an exact solution y{t) with a slightly perturbed initial condition y{0) 
(such that \y{t) — x(t)\ remain small), which stays near the numerical solution x{t) 
for sufficiently long time. We demonstrate the theorem schematically in Fig. 16.81 

To summarize, as we have mentioned in Section 16. 3[ the sensitive dependence 
of the solution of the chaotic system on the initial conditions leaves no hope to 
compute even approximate solutions with a given initial conditions. At the times 
larger than inverse of the Lyapunov exponent the numerical solution will diverge 
from the exact one. However, the shadowing theorem provides a different sense 
of numerical solutions for ODEs. It assures that numerical solutions of dynamical 
systems are meaningful. Usually one is not interested in a particular solution of 
the dynamical system in the chaotic regime but rather in the structure of attractors 
that evolve from a set of initial data. Since numerically obtained solutions are close 
to some exact solutions with close initial values, we can use them to understand 
general characteristics of a given dynamical system such as its phase portrait. 
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Figure 6.8: Shadowing theorem: x{t) is an exact solution, x{t) is numerically ob- 
tained solution that diverges from x{t) exponentially. y{t) is an exact 
solution of the same system but with slightly different initial conditions. 
y{t) "shadows" x{t). 



CHAPTER 7 

Renormalization in FPU chains: analytical description 



7.1 Hamiltonian formulation of FPU chains 

Consider a chain of particles coupled via nonlinear springs as shown in Fig. 11.11 
Suppose the total number of particles is and the momentum and displacement 
from the equilibrium position of the j-th particle are pj and qj, respectively. In this 
thesis, we consider the systems with only the nearest-neighbor interactions. Then 
the chain can be described by the Hamiltonian 

H = H2 + V. (7.1) 

The quadratic part of the Hamiltonian takes the form 

1 ^ 

H2 = ^J2p" + iq,-q,~ir, (7.2) 
and the anharmonic potential V is the function of the relative displacement {qj — 

N 

i=i 

where v{Aq) is a potential of the spring between two adjacent particles with the 
distance Aq between them. Here periodic boundary conditions g_i = qN and p-i = 
Pn are imposed. Since the total momentum of the system is conserved, it can be set 
to zero. In this thesis, we only consider the potentials of the restoring type, i.e., the 
potentials for which the Gibbs measure exists. In order to study the distribution of 
energy among the wave modes, we transform the Hamiltonian fl7.ip to the Fourier 
variables Qk, Pk via Eq. (12. 2p . In Chapter [21 it was shown that this transformation 
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is canonical. The Hamiltonian (17.11) becomes 

1 



fc=i 



where ujk = 2sm{7Tk/N) is the hnear dispersion relation. For the case of the /3-FPU 
chain, when only the fourth order potential is present, we will present the form of H 
below [Eq. (19.11) ]. Note that, throughout the thesis, for the simplicity of notation, we 
denote the periodic wave number space by the set of integers in the range [0, — 1], 
i.e., we drop the conventional factor, 2tt/N. The zeroth mode of the momentum Pq 
vanishes due to the fact that the total momentum is zero. The zeroth mode of the 
displacement Qo can also be set to zero, which is the consequence of the fact that 
the Qj in the Hamiltonian (17. ip are determined up to a constant. Therefore, we have 
the following conditions 



N-l 



J2p, = 0, (7.5) 



j=0 
N-l 

5]g, = 0. (7.6) 

Now let us demonstrate why the wave description is convenient on the example of 
noninteracting waves, i.e., when V = 0. Here, we discuss the distribution of energy 
among the wave modes. However, in the completely linear system there is no energy 
exchange among the wave modes and, therefore, the system will never reach energy 
equipartition state. Therefore, we imagine that there exist a weak nonlinearity that 
mixes energy and thus we can talk about equilibrium. But on the other hand, we 
consider this nonlinearity to be so weak that practically all the energy is in the linear 
modes. The other way of thinking about the issue of reaching thermal equilibrium 
is to imagine a heat bath, which is in thermal contact with the chain. In this case, 
the particles of the chain can exchange energy through the thermal bath. Then, 
even though the interactions between the particles of the chain are harmonic, it 
makes perfect sense to discuss thermal equilibrium of the chain. These two ways of 
describing thermal equilibrium correspond to microcanonical and canonical (Gibbs) 
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distributions, respectively. In this Chapter, we will first discuss the microcanonical 
distribution and obtain the renormalization property of the linear dispersion of the 
FPU chains. Then, in Section 17.71 we will turn our attention to the canonical 
distribution and obtain similar results about the renormalization. Then, we will 
discuss why the canonical distribution is more convenient in practical calculations. 

7.2 Free noninteracting waves 

If the nonlinear interactions are weak, and can be neglected, Eq. (17. 4p takes 

form 



In this case, it is convenient to further transform the Hamiltonian (17.71) to the 
complex normal variables defined by Eq. (12.61) . As shown in Chapter [2], this 
transformation is also canonical, i.e., the dynamical equation of motion is given by 
Eq. (12.81) . In terms of these normal variables, the Hamiltonian (17.71) takes the form 

N-l 

H=J2^k\ak\^- (7-8) 

k=l 

For the system of noninteracting waves, it is possible to obtain a standard virial 
theorem [33] in the form 




N-l 



(7.7) 



k=l 



{Kk)\v=o = {Uk)\v=o, 



(7.9) 



where 



(7.10) 



is the kinetic energy of the k-th mode and 



Uk = ^^l\Qk 



(7.11) 
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is the quadratic potential energy part of the k-th mode. Moreover, in thermal 
equilibrium we have equipartition of energy [Eq. f l3.25p ] in the form 

{Kk) = {Ki), (7.12) 
m = {Ui). (7.13) 

To show this, one can notice that the terms \Pk\^ and uj1\Qk\^ appear in the summa- 
tion in Eq. (17. 7p as independent items. Therefore, when the averages in Eqs. (17.121) 
and (17. ISp are computed using either microcanonical or canonical measure, these 
averages hare independent of the wave numbers. Let us show this for the kinetic 
energy using the Gibbs measure 




where I{Q) is part of the integral that only depends on Q. To obtain the last 
equality in Eq. (I7.14p . the change of variables P ^ z was made in order to indicate 
that the resulting expression does not depend on the wave number k. Therefore, 
the RHS of Eq. (I7.14p is independent of the wave number and so is the average 
kinetic energy of each wave mode. In a similar way, we can show independence of 
the average potential energy of the wave mode. Surprisingly, as we will see below, 
this property holds in the thermal equilibrium even when nonlinearity is present. 
As a consequence of energy equipartition of system (17. 8p . we have the following 
properties of free waves, 

{alai) = (|a,H5f, (7.15) 
(akai) = 0, (7.16) 
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where 6 is temperature given by 

e = {\Pk\'). (7.17) 

Here, we have used the equipartition theorem that was discussed in Chapter [3] 
[Eq. (13.24] . Note that Eq. (17.151) gives the classical Rayleigh- Jeans distribution for 
the power spectrum of free waves [62] 

nk = — . (7.18) 

In order to set the grounds for the case of nonlinear FPU chains, we rewrite 
Eqs. (I7.15P and (I7.16P in terms of Fourier modes. Let us invert Eq. (12. 6p 

Pk = y^(a. + a^_,), (7.19) 
Qk = y^(«. + aW). (7.20) 
Then, we can compute the following correlators using Eqs. (I7.15P and (I7.16P 



iPkPi) = {\p'M^\ 

■k+l 
N ; 



(QkQi) = {\Qk\')S''+' 



(PkQi) = 0. 

(7.21) 

Next, we substitute Eq. into Eqs. frTTSD and fl7:T6|l and then use Eq. (\T2l} to 
obtain 

(«» = 7^m\')+u;l{\Qk\'))6f = -6f, (7.22) 



2uJk oJk 
1 ^ 

2ujk 



M = -^((|Pfc|2)-a.^(|g,p))4+' = o. (7.23) 



The transformations that we just did may seem trivial since in the linear regime 
Eqs. (I7.15P and (17.161) are sufficient to give a statistical description of different 
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modes. However, situation changes when nonhnearity is present. In particular, the 
waves Ofc and a^-k become correlated, i.e.. 



since the property (17. 9p is no longer valid. However, Eq. (17.211) stays valid even in 
the nonlinear regime and the proof consequences of this fact will be discussed in the 
next Sections. 

7.3 Nonlinear interactions and microcanonical description 

Now, we turn to the case of nonhnearity of any strength. As we have mentioned 
above, the waves do not constitute a set of uncorrelated waves. However, as we 
will show below, a complete set of new renormalized variables hk can be constructed. 
Using these new variables, the strongly nonlinear system can be viewed as a system 
of "free" waves. These waves are free in the sense of vanishing correlations and the 
power spectrum, i.e., the new variables (ik satisfy the properties of free waves given in 
Eqs. (I7.15P and (17.161) . Next, we show how to construct these renormalized variables 
(ik- As we have discussed in Chapter [3l the microcanonical measure [Eq. 13. 6j (total 
energy and the number of particles are fixed), can be used in order to study the 
statistics of the system (17. ip 



where E is the total energy and dpdq = dpi . . . dp^dqi . . . dq^. Let us make a change 
of variables 



(7.24) 




(7.25) 



(gi, . . . , gAr) ^ {yi = qi,y2 = q2 - qi, ■ ■ ■ ,yN = qN - qN-i) 



(7.26) 
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This transformation is non-degenerate since it is given by the matrix 

( \ . . ^ 
-11.. 



M 



\ 



. . 1 
. . -1 1 



(7.27) 



with the unity determinant. Under the transformation ( 17.26p . the microcanonical 
measure f l7.25p becomes 



dw{y, y) = S{H{p, y) - E)dp2 ■ ■ ■ dpNdy2 ■■■dy_ 



N, 



(7.28) 



where the Hamiltonian takes the form 



1 ^ 



i=2 



^ /I 

i=2 



+ Pn) 



+ 5(y2 + 



■■ + !/«)' + i'(-(!/2 + •■■ + !/«)) (7.29) 



Note that the measure fl7.28l) does not prescribe any probabihty for pi and yi, 
which is a reflection of the fact that in the coordinates (p, q) the center of mass is 
at rest at the origin [Eqs. (17. 5p and (17. 6p ]. Therefore, we have lowered the num- 
ber of degrees of freedom by 2. Although the transformation to the new variables 
yj is non-degenerate, it is not canonical. It means that the pair {p, y) does not 
a constitute a canonically conjugate pair of variables. However, we will only use 
the non-degeneracy condition below. Next, we will use the properties of Hamilto- 
nian ( 17.29P in order to study the statistical properties of pj and yj. 



7.4 Statistics of oscillators 

Since Hamiltonian (17.290 is an even function of p, we obtain that the correlation 
between any ps and any yj vanishes as an integral of an odd function over the whole 



46 



space with respect to measure fl7.28p . i.e., 

M.-fp^yM-\P^-A-l-P,yM. (7.30) 



where {. . .)m denotes a microcanonical averaging. If the number is equal to its 
negative, than it must be zero 

Mm = 0. (7.31) 

In order to study the statistical properties of yj it is convenient to introduce another 
variable iji 

yi = -y2 Vn, (7.32) 

or, in terms of qj, 

yi = Qi~ Qn- (7.33) 

Since the Hamiltonian (17.291) is independent of yi, as we discussed above, let us 
abuse the notation and drop the "check" sign in yi. From now on, we have 

yi = yi. (7.34) 

However, it is important to emphasize again that the microcanonical measure ( 17.28^ 
only prescribes the probabilities for pj and yj with j G [2, A^] and pi and yi are 
determined from Eqs. (17.51) and (17.61) . To symmetrize the measure dw, we introduce 
the microcanonical measure over the whole set of indices j G [1, A^] via 

dwi{p,y) = 5{H{p,y) - E)5 ^S^^jJ ^ (^^^) ' ' ' ^^^^^^ ' ' ' '^'^'^^^ 
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We can view the measure dw as a projection of dwi on to the space of pj and yj 
with J G [2, N] 

dw = / dwi. (7.36) 

Using the symmetry properties of dwi, below we will obtain an important property 
of p and y, i.e., {yl)m and {ykyi)m are independent of k and I. Therefore, we define 

{y^)m = {yl)m, (7.37) 
{yy)m = {ykyi)m, k^l. (7.38) 

With these new definitions, we can rewrite the correlation 

{ykyi)m = {y^)mSi + {yy)mi^ - ^i), (7-39) 
for any k and /. Using Eqs. f l7.32p and (17.341) . we obtain for 7^ 1 

{yy)m = {yi,yk)m = {-y2 yN, yk)m = -{yy)m{,N - 2) - (7.40) 

Thus, we find the following relation between {y'^)m and {yy)m 

{yy)m = -j^{y^)m. (7.41) 

Note that in the thermodynamic limit (canonical ensemble), when oo, we have 

{yy)m = 0. (7.42) 

Similarly, we find the following equation for p 

{pp)m = -j;^{p')m- (7.43) 
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To summarize, we have obtained the following expressions for the correlations of p 
and y 

{PkPl)m = {p')mSf-^{l-Sf), (7.44) 

{ykyi)m = {y')mSf-^^ii-6f), (7.45) 

{Pkyi)m = 0. (7.46) 



We will use these properties in order to study the statistics of the Fourier modes 
and, after that, of the normal modes. 



7.5 Statistics of Fourier modes 

The Fourier transform of yj is defined by 

^^ = 7^E^^-^^' (7.47) 

with Fq = as a consequence Eq. (I7.32p . Next, let us use Eq. (I7.45P to obtain the 
statistical properties of Y^. 

/ AT N \ N 

{YkYi)m = {j;^2^yje^ l^ysC ^ ) =^2^{yjys)me ^ 
\ i=i s=i I ^ j,s=i 



N -I 

(7.48) 



for k,l ^ N . Here, we have used the following equality 

/ 27rirn\ _ ,^ ^. 

exp ( ) = ^^N. (7.49) 



n=l 
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for any integer r and j. The sum is equal to 1 if r is a multiple of A^. Otherwise, 
the sum is zero. Equation fl7.49p is just a summation of the geometrical progression 
with the base exp(— 27rir/A^). Similarly, using Eq. fl7.44p . we find for P 

{PkPi)m = j^{p')^5'i,^'. (7.50) 

Finally, from Eq. (17.461) . it immediately follows that 

{PkYi)m = 0. (7.51) 

Next, we find the relation between and using Eqs. (12.21) and (17.471) 

N N N 

i ^ — \ 2nikj i ^ — \ 27rifcj i V ■\ 2irikj 



Qk - <5fce'^' = -lexp ( ^ ) LUkQk- 



(7.52) 



Now, we have the following statistical properties of the Fourier modes Pk and Qk 

{PkPl)m = J^y)m5%^\ (7.53) 

{QkQi)m = \^^{y')m5%^\ (7.54) 

{PkQi)m = 0. (7.55) 

In the next Section, we will use these properties to study the statistics of the normal 
modes. 

7.6 Renormalized waves and their statistics 

As we have seen in Section 17.21 if the anharmonic part of the potential is 
sufficiently weak, then corresponding waves ak remain almost free, and Eqs. (17.151) 
and (I7.16P would be approximately satisfied in the weakly nonlinear regime. How- 
ever, when the nonlinearity becomes stronger, waves become strongly correlated. 
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and, in general, the correlations between waves [Eq. fl7.16p ] no longer vanish. In 
particular, {auaN-k) 7^ 0. Naturally, the question arises: can the strongly nonlinear 
system in thermal equilibrium still be viewed as a system of almost free waves in 
some statistical sense? In this thesis, we give an affirmative answer to this ques- 
tion. It turns out that the system (17. ip can be described by a complete set of 
renormalized canonical variables a^, which still possess the wave properties given 
by Eqs. (I7.15P and (I7.16P with a renormalized linear dispersion. The waves that 
correspond to these new variables hk will be referred to as renormalized waves. We 
will show that these renormalized waves possess the equilibrium Rayleigh- Jeans dis- 
tribution [62] and vanishing correlations between waves. Therefore, they resemble 
free, non-interacting waves, and can be viewed as statistical normal modes. Fur- 
thermore, it will be demonstrated that the renormalized linear dispersion for these 
renormalized waves has the form Uk = ri{k)ujk, where ri{k) is the linear frequency 
renormalization factor. Moreover and quite surprisingly, this renormalization fac- 
tor is independent of /c as a consequence of the microcanonical measure (or Gibbs 
measure in the thermodynamic limit). 

Consider the generalization of the transformation (12. 6p . namely, the trans- 
formation from the Fourier variables Qk and to the renormalized variables 
by 

Ofc = 7= — , (7.56) 

where Uk is an arbitrary positive function, which satisfies Eq. (12. 7p . In Chapter [2] 
it was shown that transformation (I7.56P is canonical. From Eqs. (I7.53p . (I7.54p . 
and (I7.55p . we obtain the following correlator of 5^ 

{0'kO'l)m = -Z—7=={{PkPl)m—'^k^l{QkQl)m) 

^ '{p'U-^iy'uV^^^. (7.57) 



2y/uJku;i N -1\ LJ^ 
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Let us choose ujk such that {akai)m vanishes for any k and /. For this, we demand 



{P')m-4(y')m = (7.58) 



Thus, we obtain 



^OJk_ I {P^)ra _ {K)rn 

^ - o;. - V {y^U ~ V {UU ' ^ ^ ^ 

where 7] is the renormahzation factor, K is the kinetic energy 

1 ^ 
i=i 

and U is the quadratic potential energy 

1 ^ 

Note that the RHS of Eq. fl7.59p is independent of the wave number k. Therefore, 
we obtain that the renormahzation factor rj is also independent of the wave number 
k. For the power spectrum, we obtain 

{ak~a:)rn = ({P\n + ^'U^f = T^^^f- (7-62) 

2^ijj\JjJi iV — 1 V LO^ / iV — 1 Wfc 

To summarize, we have the following properties of the renormalized waves 

{akai)m = 0, (7.63) 

{ak~al)m = #T^^'- ^^-64) 
TV - 1 Uk 

In our method, the construction of the renormalized variables Ofc does not de- 
pend on a particular form or strength of the anharmonic potential, as long as it is 
of the restoring type with only the nearest-neighbor interactions, as in Eq. (17. ip . 
Therefore, our approach is non-perturbative and can be applied to a large class 
of systems with strong nonlinearity. However, in this thesis, we will focus on the 
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/3-FPU chain to illustrate the theoretical framework of the renormalized waves. 
We will verify that effectively constitute normal modes for the /3-FPU chain 
in thermal equilibrium by showing that (i) the theoretically obtained renormalized 
linear dispersion relationship is in excellent agreement with its dynamical mani- 
festation in our numerical simulation, and (ii) the equilibrium distribution of 
is still a Rayleigh- Jeans distribution and a^s are uncorrelated. Note that similar 
expressions for the renormalization factor 77 have been previously discussed in the 
framework of an approximate virial theorem [2] or effective long wave dynamics via 
the Zwanzig-Mori projection [35]. However, in our theory, the exact formula for 
the renormalization factor is derived from a precise mathematical construction of 
statistical normal modes, and is valid for all wave modes k — no longer restricted 
to long waves. 

7.7 Canonical ensemble and Gibbs measure 

We have used the microcanonical measure in order to describe the statistical 
properties of the normal modes of the FPU chain. However, this approach has at 
least two disadvantages. The first one is that the microcanonical measure is hard to 
deal with, in particular, it is hard to compute {K)^ and (f/)^ and hence to obtain 
the value of the renormalization coefficient rj [Eq. (17.591) ]. And the second disad- 
vantage of the microcanonical description is that the isolated system of the coupled 
oscillators does not necessarily reach an equilibrium state as it was observed in the 
experiment by FPU [18]. The nonlinear potential V in Eq. (17.11) should be strong 
enough to "drive" the chain to chaos. Both of these drawbacks can be overcome if 
we consider the chain of oscillators to be in equilibrium with an imaginary thermal 
bath, i.e., we describe the chain as a canonical ensemble. Then the Gibbs measure 
provides a statistical description of the chain and the Gibbs measure is much eas- 
ier to compute with. Furthermore, we do not have to worry about the system not 
reaching thermal equilibrium since the interactions with the thermal bath equili- 
brate the chain. As we have noted in Chapter [3l both microcanonical and canonical 
descriptions coincide for large systems, i.e., when ^ 00. 

Consider the FPU chain in equilibrium with the thermal bath (thermodynamic 



53 



limit, N oo). Then, we can write the Hamiltonian fl7.29p as 

N ^ N 

The Gibbs measure is given by Eq. fl3.20p . with the partition function 

Z = j exp (^^H{p, y)^ dpdy. (7.66) 

Here, we integrate over the whole set of variables pj and yj for j G [l,iV]. Note 
that here we use the notation 6 for temperature since we will use the notation T 
for the interaction coefficient below. Now, the probability measure dw{p,y) can be 
written as a product of probability measures for each individual component pj and 
yj since exponential of the sum is equal to the product of the exponentials of the 
corresponding items of the sum 

N N 

dw{p, y) = Yl dwp{pj) Yl dwy{yj), (7.67) 

where 

->2 



and 



Random variables A and B are called independent if 

P{A and B) = P{A)P{B). (7.70) 



Therefore, from Eq. (17.671) . we see that all pj and all yj form a set of 2N independent 
random variables. 
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We compute the normalizing constants and Zy 



Z, = I exp[-^]dp= V2^. (7.71) 



Similarly, 



Zy = exp (^-i (^^ + v{y)^^ dy. (7.72) 

Note that 

Z = (ZpZyf. (7.73) 

From the independence of all momenta and all displacements, we obtain the follow- 
ing relationships 

{PkPi) = {p')Sf, {ykyi) = {y')Sf, {pm) = 0, (7.74) 

where (. . .) denotes averaging with respect to the Gibbs measure (13.201) . It is in- 
structive to compare Eq. fl7.74p with the similar properties given by Eqs. fl7.44p . 
( I7.45p . and ( I7.46P for the microcanonical averaging. As it is expected, in the limit 
^ oo, both ensembles give the same result. We then follow the same chain of 
transformations as we did for the microcanonical ensemble. First, we transfer to the 
Fourier modes and Qk via Eqs. (12. 2p and (I7.47p . Using Eq. (17.741) . we obtain the 
following correlation properties of the Fourier modes 



(P,P,) = (7.75) 
{QkQi) = (7-76) 
(PkQi) = 0. (7.77) 



As a next step, we transfer to the renormalized variables given by Eq. (I7.56p . We 



55 



obtain the following correlators of the renormalized waves using Eqs. (I7.75p - fl7.77p 



We choose oj^ to annihilate (af^ai) for all A; and Z and obtain the same renormalization 
condition Uk = V^k as we had for the microcanonical ensemble [Eq. fl7.59p 



However, now we can actually compute rj 



(7.80) 



CO 



V=iOZyr^\^J y'e''^'^"''^J dyj (7.81) 

Furthermore, the power spectrum (17.790 takes a simple form 

(S.S;) = ^5f. (7.82) 

Again, we compare Eq. fl7.82p with the similar expression for the microcanonical 
ensemble given by Eq. ( ]7.64p . Note that if we identify {p'^)m as temperature then, 
in the thermodynamic limit, both expressions produce the same result. 

The immediate consequence of the fact that rj is independent of k is that 
the power spectrum of the renormalized waves possesses the precise Rayleigh- Jeans 
distribution, i.e., 

Uk = ^, (7.83) 
u^k 

from Eq. fl7.79p . where fik = {\dk\'^). Combining Eqs. (12. 6p and (17.560 . we find the 
relation between the "bare" waves and the renormalized waves dk to be 

^^ = l{Vv + ^)a,+ ^-(^Vv-^) aN^k. (7.84) 
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Using Eq. fl7.84p . we obtain the following form of the power spectrum for the bare 
waves ttfc 

n, = Ul + \) —, (7.85) 

which is a modified Rayleigh- Jeans distribution due to the renormalization factor 
(H-l/r7^)/2. Naturally, if the nonlinearity becomes weak, we have rj 1. Therefore, 
all the variables and parameters with tildes reduce to the corresponding "bare" 
quantities. In particular, ojk ^k, Ofc — ^ Qfc, ""-fc — ^ ^fc- It is interesting to point 
out that, even in a strongly nonlinear regime, the "free- wave" form of the Rayleigh- 
Jeans distribution is satisfied exactly [Eq. fl7.83p ] by the renormalized waves. Thus, 
we have demonstrated that even in the presence of strong nonlinearity, the system 
in thermal equilibrium can still be viewed statistically as a system of "free" waves 
in the sense of vanishing correlations between waves and the power spectrum. 

Note that, in the derivation of the formula for the renormalization factor 
[Eq. f l7.59p ]. we only assumed the nearest-neighbor interactions, i.e., the potential 
is the function of qj — qj+i- One of the well-known examples of such a system is 
the /3-FPU chain, where only the forth order nonlinear term in V is present. In the 
remainder of the thesis, we will focus on the /3-FPU to illustrate the framework of 
the renormalized waves dk- 



CHAPTER 8 
Numerical study of the /3-FPU chain 



The Hamiltonian of the /?-FPU chain is of the form 

^ = E 2^? + 2^^^ - ^^+1)' + - (8.1) 

where /3 is a parameter that characterizes the strength of nonhnearity. The corre- 
sponding equations of motion become 



(8.2) 

Pj = - % + Qj+i) +P[{qj+i - Qjf - {qj - Qj-iY] ■ 

To investigate the dynamical manifestation of the renormahzed dispersion Uk of 
Sfc, we numerically integrate Eq. (18.21) . Since we study the thermal equilibrium 
state [1,9,20,36] of the /5-FPU chain, we have verified that the results discussed in 
the paper do not depend on details of the initial data: we have used random initial 
conditions, i.e., pj and qj were selected at random from the uniform distribution in 
the intervals (— Pmax,Pmax) and (— ^max, Q'max), respectively, with the two constraints 
that (i) the total momentum of the system is zero and (ii) the total energy of the 
system E is set to be a specified constant. The results we obtained in the thermal 
equilibrium state were independent of the initial condition. 



8.1 Parametrization of the FPU chains 

Note that the behavior of the FPU for fixed number of particles is fully 
characterized by only one parameter [49]. Let us show this for the Hamiltonian 
with = 1 and then generalize it for any N. Consider the general form of the 
Hamiltonian 

H=lip' + y') + ^f, (8.3) 
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where r > 2. Suppose 7 = PE'''^'^^^, where E is the total energy of the system fl8.3p . 
Consider another system with the same Hamihonian but different total energy Ei 
and the nonlinearity parameter (3i 

Hi = l{pl + yl) + %l. (8.4) 

If f^iEl^^'^ = 7 then every trajectory {pi,qi) of system (18. 4p can be obtained by 



scaling a corresponding trajectory of system (18. 3p by a factor ^jE^jE 

ipi,yi) = ^^ip,y)- (8.5) 

In order to prove this, we combine Eq. (18.51) with Eq. (18.41) 

r/2 



16) 



By comparing Eq. (18.31) with Eq. (18. 6p we obtain 



r/2-1 

which yields the following equality 

^ p^E[/^-\ (8.8) 

Now we generalize this argument for any number A^. In this case, in Hamiltoni- 
ans (18. 3p and (18.41) . we will have summation over each particle number j. Then, 
we apply exactly the same argument to each item in that summation as we did for 
= 1. In particular, for /?-FPU chain, we have r = 4 and, therefore, it can be 
fully characterized by one parameter j3E. In the numerical experiments, when we 
test our results for different strength of nonlinearity, we will hold the total energy 
E fixed and vary the nonlinearity parameter (3. 



59 



8.2 Details of the numerical experiments 

We use the sixth order symplectic Yoshida algorithm [58], which was briefly 
discussed in Chapter |5l In most of the numerical experiments, the system size was 
chosen to be = 128 or N = 256 and the total energy E = 100. We probed the 
/3-FPU chain with a wide range of the nonlinearity strength f3 G [10^^,10^]. The 
time step is chosen to be dt = 0.01, which ensures the conservation of the total 
system energy up to the ninth significant digit for a runtime r = 10® time units. 
This runtime was enough for the system to reach thermal equilibrium even for the 
small nonlinearity f3 = 10~^. 

8.3 Thermalization 

In order to confirm that the system has reached the thermal equilibrium 
state [46], the value of the energy localization [15] was monitored via 

where Gj is the energy of the j-th particle defined as 

Gj = ^p'^j + ^[{qj-qj+if + iqj-i-qj)'^] 

+ ^[{Qj-Q^ir + iQj-i-Qjn (8.10) 

If the energy of the system is concentrated around one site, then L{t) = 0{N). 
Whereas, if the energy is uniformly distributed along the chain, then L{t) = 0(1). 
In our simulations, in thermal equilibrium states, L{t) is fluctuating in the range 
of 1-3. In Fig. 18.11 we plot the time evolution of L{t) for = 256, E = 100, and 
(3 = 1. 

The energy localization function L{t) is also used in detecting the appearance 
of discrete breathers — spatially localized periodic solutions of the discrete lattice 
that we will discuss in Chapter [T2l 

Since our simulation is of microcanonical ensemble, we have monitored various 
statistics of the system to verify that the thermal equilibrium state that is consistent 
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Figure 8.1: Energy localization L(t) computed via Eq. (18.91) . The chain was mod- 
eled for = 128, /5 = 1, and E = 100. 



with the Gibbs distribution (canonical ensemble) has been reached. Moreover, we 
verified that, for as small as 32 and up to as large as 1024, the equilibrium 
distribution in the thermalized state in our microcanonical ensemble simulation is 
consistent with the Gibbs measure. We compared the renormalization factor (17.590 
by computing the values of {K) and (f/) numerically and theoretically using the 
Gibbs measure and found the discrepancy of 7] to be within 0.1% for (3 = 1 and the 
energy density E /N = 0.5 for from 32 to 1024. 



8.4 Spatiotemporal spectrum 

We now address numerically how the renormalized linear dispersion uj^ mani- 
fests itself in the dynamics of the /3-FPU system. We compute the spatiotemporal 
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Figure 8.2: The spatiotemporal spectrum |dfc(c<j)p in thermal equihbrium. The chain 
was modeled for = 256, j3 = 0.5, and E = 100. [max{— 8, In \ak{uj)\'^} , 
with corresponding gray scale, is plotted for a clear presentation]. The 
solid curve corresponds to the usual linear dispersion uj^ = 2 sin(7rA;/A^). 
The dashed curve shows the locations of the actual frequency peaks of 
|afc(^)p. 

spectrum |aj!c(u;)p, where ak{u) is the Fourier transform of afc(t). (Note that, for 
simplicity of notation, we drop a tilde in a^.) Figure [8l2l displays the spatiotemporal 
spectrum of dk, obtained from the simulation of the /?-FPU chain for = 256, 
f3 = 0.5, and E = 100. In order to measure the value of t] from the spatiotemporal 
spectrum, we use the following procedure. For the fixed wave number k, the corre- 
sponding renormalization factor r]{k) is determined by the location of the center of 
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the frequency spectrum \ak{uj)\'^, i.e., 

Vik) = , with ujc{k) = ... . , ■ 

The renormahzation factor r]{k) of each wave mode k is shown in Fig. 18.31 The 
numerical approximation fj to the value of rj is obtained by averaging all ri{k), i.e., 



k=l 

The renormahzation factor for the case shown in Fig. 18.21 is measured to be 77 ~ 
1.1824. It can be clearly seen in Fig. 18.31 that ri{k) is nearly independent of k and its 
variations around f] are less then 0.3%. We also compare the renormahzation factor 
T] obtained from Eq. (I7.59P (solid line in Fig. 18.31) with its numerically computed 
approximation f/ (dashed line in Fig. 18.31) . Equation (17.591) gives the value r] ^ 1.1812 
and the difference between rj and f] is less then 0.1%, which can be attributed to 



the statistical errors in the numerical measurement. In Fig. 18. 4[ we plot the value 
of ?7 as a function of f3 for the system with N = 256 particles and the total energy 
E = 100. The solid curve was obtained using Eq. (17.801) while the circles correspond 
to the value of t] determined via the numerical spectrum 10,^(07)^ as discussed above. 
It can be observed that there is excellent agreement between the theoretic prediction 
and numerically measured values for a wide range of the nonlinearity strength (3. 
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Figure 8.3: Independence of k of the renormalization factor i]{k). The circles cor- 
respond to ri{k) obtained from the spatiotemporal spectrum shown in 
Fig. 18.21 [only even values of k are shown for clarity of presentation] . The 
dashed line corresponds to the mean value f/. For (3 = 0.5, the mean 
value of the renormalization factor is found to be f/ ^ 1.1824. The varia- 
tions of rjk around f] are less then 0.3%. [Note the scale of the ordinate.] 
The solid line corresponds to the renormalization factor rj obtained from 
Eq. fl7.59p . For the given parameters r] ^ 1.1812. 
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Figure 8.4: The renormalization factor as a function of the nonhnearity strength 
(3. The analytical prediction [Eq. f l7.80p ] is depicted with a solid line 
and the numerical measurement is shown with circles. The chain was 
modeled for N = 256, and E = 100. 



CHAPTER 9 
Dispersion relation and resonances 

In this chapter, we will discuss how the renormalization of the linear dispersion of the 
/3-FPU chain in thermal equilibrium can be explained from the wave resonance point 
of view. This formalism is used in the theory of wave turbulence, that we briefly 
discussed in Chapter HI In order to give a wave description of the /3-FPU chain, we 
rewrite its Hamiltonian in terms of the renormalized variables dk- Then, we study the 
interactions of the renormalized waves. In order to address how the renormalized 
dispersion arises from wave interactions, we study the resonance structure of our 
nonlinear waves. We will demonstarte that the /3-FPU system is a Hamiltonian 
system with four-wave interactions. We will discuss the properties of the resonance 
manifold associated with the /3-FPU system as a first step towards the understanding 
of its long time statistical behavior. We comment that the resonance structure is 
one of the main objects of investigation in wave turbulence theory [7,39,42,43,52, 
53,62]. The theory of wave turbulence focuses on the specific type of interactions, 
namely resonant interactions, which dominate long time statistical properties of the 
system. On the other hand, the non-resonant interactions are usually shown to have 
a total vanishing average contribution to a long time dynamics. 

In order to rewrite the Hamiltonian of the /3-FPU chain given by Eq. (18.10 in 
terms of renormalized waves, we first transfer to the Fourier space using Eq. (12. 2p . 
The Hamiltonian becomes 

1 

k=l 

+ ^ E ^ku;i^n.u;s[-{QkQiQrnQ:s'J'^ + cx.) + QkQiQ*^Q:6tl 

k,l,m,s=l 

(9.1) 

Here, we have used Eq. (17.491) . Next, we transfer to the renormalized waves ak using 
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the inverse of Eq. fl7.56p 



Qk = —7={ak - a*^-k)- 

y/2uJk 



In order to give a wave description of the /3-FPU chain, we rewrite Hamiltonian (18 .ip 
in terms of the renormahzed variables dk [Eq. ( I7.56p ] with = rjUk, 



N-l 

+ E 



k=l 
N-l 

nkl 

- ms 

k,Lm,s=l 



A^^dkdid*^dl + { '^Af^'dkdidmdl + c.c. 



(9.3) 



where c.c. stands for complex conjugate, and 

3/5 

= ^^V^kUJliOn^UJs (9.4) 

is the interaction tensor coefficient. Note that, due to the discrete nature of the 
system of finite size, the wave space is periodic and, therefore, the "momentum" 
conservation is guaranteed by the following "periodic" Kronecker delta functions 

\kl _ srkl TklN /q r\ 

A Mm s.klm s.klm < s,klm (q, r;\ 

'-^s — "s "sAf "T ^sNN^ 

A fcZms _ s,klms s,klms s,klms (q 

^0 — ^NN "N "NNN- W-'J 

Here, the Kronecker 5-function is equal to 1, if the sum of all superscripts is equal 
to the sum of all subscripts, and 0, otherwise. The periodic delta functions arise 
when we make the change of indices of the form N — k ^ k in the summation in 
order to obtain Eq. (19. 3p . 

In analogy with quantum mechanics, where and a are creation and annihi- 
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Figure 9.1: Interaction process of type (2 2) given by the term ala^amds^ms- 

lation operators, we can view dl as the outgoing wave with frequency Uk and dk as 
the incoming wave with frequency Ok- Then, the nonhnear term 

dld^d^dsA'^, (9.8) 

in system (19.31) is schematically shown in Fig. 19. H and can be interpreted as the inter- 
action process of the type (2-^2), namely, two outgoing waves with wave numbers 
k and / are "created" as a result of interaction of the two incoming waves with wave 
numbers m and s. Similarly, 

^ims system (19.31) describes the interaction 
process of the type (3 —>■ 1), that is, one outgoing wave with wave number k is "cre- 
ated" as a result of interaction of the three incoming waves with wave numbers /, m, 
and s, respectively. Finally, afcaiamOsAQ^™'' describes the interaction process of the 
type (4 — > 0), i.e., all four incoming waves interact and annihilate themselves. Fur- 
thermore, the complex conjugate terms afca;*aj^a*Af^j, and a^a;*aj!^a*A°;^^ describe 
the interaction processes of the type (1 — * 3) and (0 — > 4), respectively. 

Instead of the processes with the "momentum" conservation given via the usual 
^ms^ 5^^"*, or ^q'™* functions for an infinite discrete system, the resonant processes 
of the /?-FPU chain of a finite size are constrained to the manifold given by Aj^^, 
Af"^, or Aq'™*, respectively. Next, we describe these resonant manifolds in detail. 
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As will be pointed out in Chapter [TT], there is a consequence of this finite size effect 
to the properties of the renormalized waves. 



9.1 Resonance manifolds 

The resonance manifold that corresponds to the (2 2) resonant processes in 
the discrete periodic system, therefore, is described by 



N 



k + I = m + s, 



(9.9) 



N 



where we have introduced the notation g = h, which means that 



9 = h, 

g = h + N, 

g = h-N, 



for any g and h. The first equation in system (19. 9p is the "momentum" conservation 
condition in the periodic wave number space. This "momentum" conservation comes 
from 



|A 



kl 



1. 



(9.10) 



klm I 



Note that |A^^| can assume only the value of 1 or 0. Similarly, from |A 
and IAq'™^! = 1, the resonance manifolds corresponding to the resonant processes 
of types (3 1) and (4 — > 0) are given by 



k + I + m = s, 

'^fc + '^Z + = , 



(9.11) 
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and 



k + l + m + s = 0, 

Ok + 1^1 + l^m + l^s = 0, 



(9.12) 



respectively. For the processes of type (3^1), the notation g = h means that 



9 = h, 

g = h + N, 

g = h + 2N. 



For the (4 — > 0) processes, g = h means that 



g = h + N, 
g = h + 2N, 
g = h + 3N. 



9.2 Trivial resonances of the type (2 2) 

To solve system (19. 9p . we rewrite it in a continuous form with 



X = 


k/N, 


y = 


l/N, 


z = 


m/N, 


V = 


s/N, 



(9.13) 



which are real numbers in the interval (0, 1). By recalling that Uk = 2r]sin{iTk/N), 
we have 



X + y = z + V, 

sin('7rx) + sin^Tcy) = sin(7rz) + sin(7rf ). 



(9.14) 
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Thus, any rational quartet that satisfies Eq. fl9.14p yields a solution for Eq. fl9.9p . 
There are two distinct types of solutions for Eq. (19.141) . The first one is given by 

X + y = z + V, 

(9.15) 

sin(7rx) + sin(7ry) = sin(7rz) + sin(7rf ). 

Here we show that the system (19.151) has only trivial solutions. We express v from 
the first equation in (I9.14p and insert it into the second equation 

sin(7rx) + smijiy) = sin(7rz) + sin(7r(a; + y — z)). (9.16) 

Using Prosthaphaeresis formulas, we can rewrite Eq. (19.160 as 

sin I TT — - — sin {7t{x — z)) sin {^{y — z)) = 0. (9-17) 



By recalling that < x,y, z < 1, we obtain that the only solutions for the sys- 
tem (19.151) are given by 

X = Z, \ X = 

or <^ (9.18) 
y = v, \^y = z, 

i.e., these are trivial resonances, as we mentioned above. 

9.3 Non-trivial resonances of the type (2 2) 

The second type of the resonance manifold of the (2 2)-type interaction 
processes corresponds to 

x + y = z + v±l, 

(9.19) 

sin(7rx) + sm{ny) = sin(7rz) + sin(7rt>). 

In order to solve Eq. (19.191) . we pull out v from the first equation and substitute it 
to the second equation. Then we use trigonometric addition formulas to obtain the 
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following two branches 



X + y 1 



H — arcsin(y4) + 2j, 



(9.20) 



Zl = 



2 TT 

x + y _ ^ 
2 



1 



arcsin(y4) + 2j, 



(9.21) 



Z2 = 
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where A = tan (7r(x + y)/2) cos (7r(a; — y)/2) and j is an integer. 

The second type of resonances arises from the discreteness of our model of a 
finite length, leading to non-trivial resonances. For our linear dispersion here, non- 
trivial resonances are only those resonances that involve wave numbers crossing the 
first Brillouin zone. This process is known as the Umklapp scattering in the setting 
of phonon scattering [65]. Note that, in the previous studies [6, 31] of the FPU chain 
from the wave turbulence point of view, the effects arising from the finite nature of 
the chain were not taken into account, i.e., only the limiting case of ^ oo, where 
is the system size, was considered. 

9.4 Numerical observations of the resonances 

In Fig. 19.21 we plot the solution of Eq. fl9.14p for x = k/N with the wave num- 
ber k = 90 for the system with = 256 particles (the values of k and are chosen 
merely for the purpose of illustration). We stress that all the solutions of the sys- 
tem (EmD are given by the Eqs. (ETIHl), (19:201) . and (I92ID, and that the non-trivial 
solutions arise only as a consequence of discreteness of the finite chain. The curves in 
Fig. I9.2l represent the loci of {z, y), parameterized by the fourth wave number v, i.e., 
X, y, z, and v form a resonant quartet, where z = m/N, and y = l/N. Note that the 
fourth wave number v is specified by the "momentum" conservation, i.e., the first 
equation in Eq. (I9.14p . The two straight lines in Fig. 19.21 correspond to the trivial 
solutions, as given by Eq. fl9.18p . The two curves (dotted and dashed) depict the 
non-trivial resonances. Note that the dotted part of non-trivial resonance curves 
corresponds to branch (19.201) . and the dashed part corresponds to branch (I9.2ip . 
respectively. An immediate question arises: how do these resonant structures man- 
ifest themselves in the FPU dynamics in the thermal equilibrium? By examining 
the Hamiltonian (19. 3p . we notice that the resonance will control the contribution 
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of terms like A,^^ in the long time limit. Therefore, we address the effect 

of resonance by computing long time average, i.e., (a^aj*am.as)A^^, and comparing 
this average (Fig. 19. 3p with Fig. 19. 2[ To obtain Fig. 19.31 the /5-FPU system was 
simulated with the following parameters: = 256, (3 = 0.5, E = 100, and the aver- 
aging time window r = 400ti, where ti is the longest linear period, i.e., ti = 2tt/uji. 
In Fig. 19. 3[ mode k was fixed with A; = 90 and the mode s, a function of k, I, and 
m, is obtained from the constraint k + I = m + s, i.e., l^^sl ~ ^- Note that we 
do not impose here the condition ojk + oJi = ujm + '^s, therefore, \{dla^ (imds) A'^J is 
a function of I and m. The two dark lines show the locations, where s = and, 
therefore, dld^amdsA^^ = 0. By comparing Figs. 19.21 and 19.31 it can be observed 
that the locations of the peaks of the long time average \{dld^dmds)A^^\ coincide 
with the loci of the (2 — > 2)-type resonances. This observation demonstrates that, 
indeed, there are nontrivial (2 2)-type resonances in the finite /3-FPU chain in 
thermal equilibrium. Furthermore, it can be observed in Fig. 19.31 that, in addition 
to the fact that the resonances manifest themselves as the locations of the peaks of 
)A^^|, the structure of near-resonances is reflected in the finite width of 
the peaks around the loci of the exact resonances. Note that, due to the discrete na- 
ture of the finite /3-FPU system, only those solutions x, z, and v of Eq. (19.141) . for 
which Nx^ Ny, Nz, and Nv are integers, yield solutions k, I, m, and s for Eq. (19. 9p . 
In general, the rigorous treatment of the exact integer solutions of Eq. (19.91) is not 
straightforward. For example, for = 256, we have the following two exact quartets 
k = {k, I, m, s}: 

k = {k, N/2 - k, N/2 + k,N- k}, 
k = {k, N/2 -k,N-k, N/2 + k}, 

for k < N/2, and 

k = {k, 3N/2 -k,k- N/2, N - k}, 
k = {k, 3N/2 -k,N-k,k- N/2}, 

for k > N/2. We have verified numerically that for = 256 there are no other exact 
integer solutions of Eq. (19.91) . In the analysis of the resonance width in Chapter [TT| 
we will use the fact that the number of exact non-trivial resonances [Eq. (19. 9p ] is 
significantly smaller than the total number of modes. 
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Figure 9.2: The solutions of Eq. (19.141) . The sohd straight hues correspond to the 
trivial resonances [solutions of Eq. f l9.18p ]. The solutions are shown for 
fixed X = k/N, k = 90, = 256 as the fourth wave number v scans 
from to {N — 1)/N in the resonant quartet Eq. (19.141) . The non- 
trivial resonances are described by the dotted or dashed curves. The 
dotted branch of the curves corresponds to the non-trivial resonances 
described by Eq. (I9.20p and the dashed branch corresponds to the non- 
trivial resonances described by Eq. (19.211) . 
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Figure 9.3: The long time average \ {a\a^arnO's) ^ms\ of the /5-FPU system in thermal 
equilibrium. The parameters for the FPU chain are = 256, (3 = 0.5, 
and E = 100. {alaidma,s)^ms "^^s computed for fixed k = 90. The 
darker grayscale corresponds to the larger value of (a^a^amas)A^^. The 
exact solutions of Eq. (I9.14p . which are shown in Fig. 19.21 coincide with 
the locations of the peaks of \{ald^dmas)A'^g\. Therefore, the lighter 
areas represent the near-resonance structure of the finite /3-FPU chain. 
(The two dark lines show the locations, where s = and, therefore, 
dld^dradgA'^^ = 0.) [max{2, ln(|(a^a;*ama5)A^^|)} with the correspond- 
ing color-scale is plotted for a clean presentation]. 

9.5 Near-resonances 

The broadening of the resonance peaks in Fig. I9.3l suggests that, to capture the 
near-resonances for characterizing long time statistical behavior of the /5-FPU sys- 
tem in thermal equilibrium, instead of Eq. (19.91) . one needs to consider the following 
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effective system 



k + I = m + s, 



(9.22) 



UJk + UJi - iOm - OJsl < At^, 



where < Aco <^ uok for any A;, and Atu characterizes the resonance width, which 
results from the near-resonance structure. Clearly, Acj is related to the broadening 
of the spectral peak of each wave ^^(t) with a = k,l, m, or s in the quartet, and this 
broadening effect will be studied in detail in Chapter [HI Note that the structure 
of near-resonances is a common characteristic of many periodic discrete nonlinear 
wave systems [14,40,50]. 

9.6 Resonances of the type (3 1) and (4 0) 

Next, we prove that there are no exact (3 l)-type resonances [Eq. (19. lip 
in the /3-FPU chain. We change variables according to Eq. (19.131) to obtain 



with < {x, y, z, v} < 1. The first equation in (19.231) implies that either x+y+z = v, 
OTx + y + z = v + l,oTx + y + z = v + 2. Since the treatment of all three cases is 
similar, we only consider the second one in details 



Denote D = \ sin(7r(x + y))\. Then, using the trigonometric addition formulas, the 
properties of the modulus, and the fact that sin(7rx) > for a; G (0, 1), we obtain 




sin('7rx) + sin(7r?/) + sin('7r2;) = sin(7rf ) 



(9.23) 




(9.24) 
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for the right-hand side of Eq. (19.241) 

D = I sin(7r(a; + ?/))! 

< I sin(7rx) 1 1 cos(7r?/) | + | sin(7r?/) 1 1 cos(7rx) | 

< sin(7rx) + sin(7ry). (9.25) 

Now, consider the left-hand side of Eq. (I9.24p 

sin(7rf ) — sin(7r2;) = sin(7r(x + y + z — 1)) ~ sin(7r2;) 
= sin(7r(x + y — 1)) cos(7r2;) 

-|- sin(7r2) cos(7r(x + y — 1)) — sin(7r2;) 
< I sin(7r(x + y ~ 1)) 1 1 cos(7r2) | 

-|- I sin(7rz)|| cos(7r(x + y)) \ — sin(7r2;) 
<\sm{7c{x + y))\ = D, (9.26) 

where the use is made of | sin(a — 7r)| = | sin(a)|. Combining Eqs. (19.251) and (I9.26p . 
we obtain 

sin(7rt') — sin(7r2:) < D < sin(7ra;) + sin(7r?/). (9.27) 

From inequahty (I9.27p . it follows that Eq. (I9.24p has no solutions, and, therefore, 
there are no exact resonances of type (3 — * 1) in the /3-FPU chain. Now it is 
apparent that all the nonlinear terms dldidmas^fms non-resonant and their long 
time average (a^a«am.as)Af^^ vanishes. As for the resonances of type (4 0), 
since the dispersion relation is non-negative, one can immediately conclude that the 
solution of the system (19.120 consists only of zero modes. Therefore, the processes 
of type (4 — s> 0) are also non-resonant, giving rise to {dkdidmds) ^o""^^ = 0. In this 
thesis, we will neglect the higher order effects of the near-resonances of the types 
(3 ^ 1) and (4 ^ 0). 

In the following Chapters, we will study the effects of the resonant terms of 
type (2^2), namely, the linear dispersion renormalization and the broadening of 
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the frequency peaks of afc(t). It turns out that, the former is related to the trivial 
resonance of type (2 — > 2) and the latter is related to the near-resonances, as will 
be seen below. 



CHAPTER 10 
Self-consistency approach to frequency renormalization 



We now turn to the discussion of how the trivial resonances give rise to the dispersion 
renormahzation. This question was examined in [22] before. There, it was shown 
that the renormahzation of the hnear dispersion of the /5-FPU chain arises due to the 
collective effect of the nonlinearity. In particular, the trivial resonant interactions 
of type (2^2), i.e., the solutions of Eq. flQ.lSp . enhance the linear dispersion (the 
renormalized dispersion relation takes the form uj^ = rjUk with rj > 1), and effectively 
weakens the nonlinear interactions. Here, we further address this issue and present a 
self-consistency argument to arrive at an approximation rjsc for the renormalization 
factor 7]. As will be seen below, the self-consistency argument essentially is of a 
mean- field type, i.e., the renormalization arises from the scattering of a wave by a 
mean-background of waves in thermal equilibrium via trivial resonant interactions. 
We note that our self-consistency, mean-field argument is not limited to the weak 
nonlinearity. Very good agreement of the renormalization factor 7] and its dynamical 
approximation r^^c — for weakly as well as strongly nonlinear waves — confirms that 
the renormalization is, indeed, a direct consequence of the trivial resonances. 



10.1 Mean-field approximation of the linear dispersion 

As it was mentioned above, the contribution of the non-resonant terms have 
a vanishing long time effect to the statistical properties of the system, therefore, in 
our self-consistent approach, we ignore these non-resonant terms. By removing the 
non-resonant terms and using the canonical transformation 

ak = ^ — , (10.1) 
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where rjsc is a factor to be determined, we arrive at a simplified effective Hamiltonian 
from Eq. (19.31) for the finite /5-FPU system 



N-l , ^ V 7V-1 

^eff = 5ZY(r?sc+— 7^^iA™s«X«m«s. (10.2) 

k=\ ^ k,l,m,s=l 

The "off-diagonal" quadratic terms afcOiv-fc from Eq. (19. 3p are not present in Eq. (110.21) . 
since are chosen so that {cikaN^k) = (see Chapter [3, Eqs. (I7.15P and (I7.16P ). 
Note that Hamiltonian (110. 2p is written in a standard four-wave interaction form 
that is widely used in wave turbulence [39]. However, for the discrete system like 
FPU one should take into account Umklapp processes that we discussed in Sec- 
tion [921 

Next we will obtain an approximation to the frequency renormalizing factor r] 
as a consequence of the trivial resonance interactions. As we have discussed above, 
the nonlinear terms in Eq. (1 10. 20 

d/^CLi dfnO's^rnsi 

can be interpreted as interaction processes when two waves and ai are created as 
a result of the interaction of two other waves and [Fig. 19.1] . From Eq. (I9.18p . 
we obtain that the trivial resonant interactions are characterized by the conditions 



:io.3) 



In Fig. 110.11 we show schematically the trivial resonant interaction process. The 
contribution of the trivial resonances in iJefr is 

N-l 

^r = 4 5^r,^/|a,na,p, (10.4) 

k,l=l 

which can be "linearized" in the sense that averaging the coefficient in front of \ak\^ 




k =m 
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Figure 10.1: Trivial resonant interactions given by Eq. (110. 3p . 
in H^^ gives rise to a quadratic form 



N-l / N-l 



k=i \ 1=1 



Note that the subscript 2 in H2 emphasizes the fact that H2 now can be viewed as 
a Hamihonian for the free waves with the famihar effective hnear dispersion [22, 62] 



N-l 



^k = ^Y.T',R\~ai?). (10.5) 
1=1 

This hnearization is essentially a mean-field approximation, since the long-time av- 
erage of trivial resonances in Eq. (110.41) is approximated by the interaction of waves 
Ofc with background waves (|a/|). The self consistency condition, which determines 
rjsci can be imposed as follows. The quadratic part of the Hamiltonian (110. 2p . com- 
bined with the "hnearized" quadratic part, H^, of the quartic H^^ , should be equal 
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to an effective quadratic Hamiltonian for the renormalized waves, i.e., 

N-1 

^2 = J^cufclfifcll (10.6) 



k=l 



Therefore, the equation for the renormahzation parameter rjsc via the self consistency 
argument becomes 



AT-l y N N-1 / N-1 \ N-1 

k=i ^ k=i \ 1=1 / k=i 



~ I ~ |2 



(10.7) 



where Uk is the renormahzed hnear dispersion, which is used in the definition of our 
renormalized wave, Eq. (110.11) . and Uk = Vsc^k- Equating the coefficients of cufclcifep 
on both sides of Eq. (110.71) for every wave number k yields 

^ f M . 3/5 ,2, 

where use is made of Eq. (19.41) . After algebraic simphfication, we have the following 
equation for rjsc 

^L-^.c=f $^^/(is,r). (10.8) 

1=1 

Using the property (17.791) of the renormalized normal variables a^, we find the 
following dependence of (IfifcP) on rjsc, 

= ^ {{\m + v>fW)) ■ (10.9) 

Combining Eqs. ffTirSl) and iHJiM leads to 



r]t, - Aril -B = 0, 



(10.10) 
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where 



1=1 

N-l 



1=1 

Here, (U) and (K) are the averaged quadratic potential and kinetic energies of the 
system. The only physically relevant solution of Eq. (110. lOp is 



Vsc=\ ^ • (10.11) 



The constants A and B can be easily derived using the Gibbs measure 

A = ^ + ^{y') = ^ + — - — A — , (10.12) 

2 2/exp(-^(y2 + /?3^)j dy 

B = f(P^) = ^- (10.13) 
See Section 17.71 for more details about computations with the Gibbs measure. 



10.2 Limiting behavior of t] 

We have studied the frequency renormalization factor r] and its approximation 
via the self-consistency argument rjsc- Here, we compare the behavior of both rj 
and r]sc in the weakly nonlinear (small (3) and strongly nonlinear (large (3) limits. 
In order to study the limiting behavior of rj and r/sc, we use the canonical Gibbs 
measure as discussed in Section 17.71 There, we obtained the following expressions 
for the pdf's for the momentum pj and displacement Uj. Any pj is distributed with 
the Gaussian pdf 

dwp = Zpexp (-^^ , (10.14) 
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and any yj is distributed with the pdf 

dWy = Zye^V + ^t)) ' 

where Zp and Zy are the normahzing constants. The renormahzation factor r] of the 
/3-FPU system in thermal equihbrium is given by Eq. f l7.59p . and its approximation 
via the self-consistency argument 7]^^ is given by Eq. (110. lip . We will use the follow- 
ing expressions for the average density of kinetic, quadratic potential and quartic 
potential energies per degree of freedom of the system, defined in Eqs. (17. Sp . (I7.60p 
and fl7:6T|l 

N 2 2 ' ^ ' 

W _ if M>=i(,.). (10.17) 

In a canonical ensemble, the temperature of a system is given by the temperature 
of the heat bath. By identifying the average energy density of the system with 
e = E/N in our simulation (a microcanonical ensemble), we can determine 6* as a 
function of e and (3 by the following equation 

l^{{K) + {U) + {V))=e. (10.19) 

We start with the case of small nonlinearity /5 — > 0. Assume that in the first order 
of the small parameter (3 the temperature has the following form 



(10.20) 



84 



where 9q = 0(1) and 6i = 0(1). We find the values of 6q and 6i using the con- 
straint (110.191) . We use the following expansions in the small parameter j3 

(10.21) 



y'^ exp 



dy 



%{4eo + {6ei-i5el)f3) + o{(3^ 



:i0.22l 



[y + 772/ ) exp 



{^00 + m - 96',)^) + 0(/52). 



(10.23) 



Then, in the first order in jS, Eq. (110. 19p becomes 



^0 + 



fv/g^(4go+(6gi-9gg)/?) 
fv^(4+(|f-3^o)/?' 



2e, 



and we obtain 6q = e and = (3/4)e^. Therefore, for the average kinetic energy 
density, we have 



(10.24) 



and, for the average quadratic potential energy density, we have 



{U) l^Ve'^{Ad^ + {<6d,-l^dl)(3) I 9 



N 



e- -e'[3 + 0{p'). (10.25) 



Finally, we obtain that for small j3 



(10.26) 
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Similarly, from Eq. fllO.lip . we find the small f3 limit of the approximation 77^^ 

Vsc = l + le(3 + 0{(3^). (10.27) 

Now, we consider the case of strong nonlinearity f3 00. From Eq. (110. 19p . we con- 
clude that temperature in the large /3 hmit, which we denote as 6^00, stays bounded, 
i.e., < 6^00 < 2e, and, in the limit of large /3, we obtain for Eq. (110.191) 

^00 + 7^ ^ = 2e. (10.28) 



After performing the integration, we obtain 6^0 = (4/3)e, and the average kinetic 
energy density becomes {K)/N = (2/3)e. For the average quadratic potential energy 
density, we have 

ju) _ 1 JZy^-^^p j-ity') iy _ r(|) Uey 

For the renormalization factor, we obtain the following large (3 scaling 



Similarly, for the approximation of rjsc-, we obtain A = C^fW^^ B = 4e/3, and C = 
2V3r(3/4)/r(l/4). Therefore, the large p scaling of rjsc becomes 



V^c=\ ^ e4/34, (10.31) 



which yields Eq. (HJim . 

10.3 Comparison of 77 and rjsc 

Next, we compare the renormalization factor rj [Eq. (17.591) ] with its approxi- 
mation rjsc [Eq. (110.111) ] from the self-consistency argument. We have shown that for 
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the case of small nonlinearity, both rj and rise have the same asymptotic behavior in 
the first order of the small parameter j3 given by Eqs. fll0.26p and fll0.27p . Moreover, 
in the case of strong nonlinearity (3 —>■ oo, both rj and rise are given by Eqs. (110.301) 
and ffTO^TD . i.e., 

V-^Vsc-^ (10.32) 

Note that, in [22], we numerically obtained the scaling r] ~ which differs 

from the exact analytical result (110. 32p due to statistical errors in the numerical 
estimate of the power. In Fig. 110. 2[ we plot the renormalization factor r] and its 
approximation r]sc for the case of small nonlinearity (3 for the system with = 256 
particles and total energy E = 100. The solid line shows r] computed via Eq. (I7.59p . 
the diamonds with the dashed line represent the approximation via Eq. (110. lip , 
and the solid circles with the dotted line correspond to the small-/5 limit (110.260 
and (110. 27p . In Fig. 110.31 we plot the renormalization factor rj and its approximation 
rjsc for the case of large nonlinearity (3 for the system with = 256 particles 
and total energy E = 100. The solid line shows rj computed via Eq. (17.590 . the 
diamonds with the dashed line represent the approximation via Eq. (110. IIP , and the 
dashed-dotted line correspond to the large-/? scaling (110. 32p . Figs. 110.21 and 110.31 
show good agreement between the renormalization factor rj and its approximation 
rjsc from the self consistency argument for a wide range of nonlinearity, from (3 ~ 
10~^ to /3 ~ 10^. This agreement demonstrates, that (i) the effect of the linear 
dispersion renormalization, indeed, arises mainly from the trivial four-wave resonant 
interactions, and (ii) our self-consistency, mean-field argument is not restricted to 
small nonlinearity. 

10.4 Effective nonlinearity 

Here, we study the decrease of the effective nonlinearity of the system due to 
the dispersion relation renormalization. To measure the nonlinearity, we compute 
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Figure 10.2: The renormalization factor as a function of the nonhnearity strength 
(3 for small values of (3. The renormalization factor r] [Eq. f l7.59p ] is 
shown with the solid line. The approximation r/^c [Eq- fll0.1ip ](via the 
self-consistency argument) is depicted with diamonds connected with 
the dashed line. The small-/3 limit [Eq. (110.261) ] is shown with the 



solid circles connected with the dotted line, 
logarithmic scale. 



Note that, abscissa is of 



the following ratio 



However, when nonhnearity is strong, 33 does not provide a physically meaningful 
measure for the strength of nonlinear interactions since the trivial resonant interac- 
tions contribute towards the renormalized linear dynamics, described by Eq. (110. 6p . 
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Figure 10.3: The renormalization factor as a function of the nonhnearity strength 
(3 for large values of (3. The renormalization factor t] [Eq. f l7.59p ] is 
shown with the solid line, r^^c [Eq. fllO.llI) ] is depicted with diamonds 
connected with the dashed line. The large-/3 scaling [Eq. (110. 32p ] is 
shown with the dashed-dotted line. Note that, the plot is of log- log 
scale. 



Therefore, we introduce the renormalized measure of effective nonhnearity 



(10.34) 



where = E — H2. Using the definition of the frequency renormalization factor rj 
[Eq. (17.591) ]. we obtain 



(10.35) 
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Figure 10.4: Nonlinearity strength (black circles) and effective nonlinearity 
strength ae (red diamonds) as functions of jS. The number of parti- 
cles is = 256, the total energy is E = 100. Note that, the plot is of 
log-linear scale. 

Therefore, we can compute se and ae using Eqs. (110. 16p and (110. 17p . In Fig. 110. 4[ 
we present the dependance of both ae and ae on the nonlinearity parameter /3. The 
total energy was held fixed at the value E = 100 and (3 was varied. We observe that 
as (3 grows, se reaches the value of ~ 0.5 whereas ae only increases up to ~ 0.25. 
It demonstrates that the effective nonlinearity is still quite small even when the 
nonlinearity parameter (3 and the bare nonlinearity ae are large. 



CHAPTER 11 
Resonance width 



We further study the properties of these renormahzed waves by investigating how 
long these waves are coherent, i.e., what their frequency widths are. Therefore, we 
consider near-resonant interactions [27,38] of the renormahzed waves a^. These 
are the interactions that occur in the vicinity of the resonance manifold. We con- 
sider these near-resonances since most of the exact resonant interactions are trivial, 
i.e., with no momentum exchanges, and they, cannot effectively redistribute energy 
among the wave modes. 

We will demonstrate that near-resonant interactions of the renormahzed waves 
dk provide a mechanism for effective energy exchanges among different wave modes. 
Taking into account the near-resonant interactions, we will study analytically the 
frequency peak broadening of the renormahzed waves dk by employing a multiple 
time-scale, statistical averaging method. Here, we will arrive at a theoretical predic- 
tion of the spatiotemporal spectrum |ajt(u;)p, where dk{uj) is the Fourier transform of 
the normal variable dk{t), and uj is the frequency. The predicted width of frequency 
peaks is found to be in good agreement with its numerically measured values. 

In addition, for a finite /3-FPU chain, we will mention the consequence of the 
Umklapp scattering (see Section [9l2|) . to the correlation times of waves. 

11.1 Effective Hamiltonian 

In the Hamiltonian (110.21) . the nonlinear terms corresponding to the trivial 
resonances have been absorbed into the quadratic part via the effective renormahzed 
dispersion u^- Therefore, the new effective Hamiltonian is 

N-l N-l 

H=J2^k\ak\^+ Yl Tl^is^'^sK^iamds, (11.1) 

k=l k,l,m,s=l 
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where 



3/3 

Tms = = OAT o V'^kUJl'^rni^s, k ^ 171, and ky^S 



(11.2) 



^ms = 0) otherwise. 

The new interaction coefficient ensures that the terms that correspond to the in- 
teractions with trivial resonances are not doubly counted in the Hamiltonian flll.ip . 
This new interactions in the quartic terms include the exact non-trivial resonant and 
non-trivial near-resonant as well as non-resonant interactions of the (2 2)-type. 

We change the variables to the interaction picture by defining the correspond- 
ing variables bk via [39,40] 



Then, the dynamics governed by the Hamiltonian (111.11) takes the familiar form 

N-l 

4 = 2 J2 f^s^lshlhmhse'''-\ (11.3) 

i,m,s=l 

where 

'^ms = '^fc + ^/ - - t^s- (11-4) 



Without loss of generality, we consider only the case oi k < N/2. As we have 
noted before, only for a very small number of quartets does oj^^ vanish exactly, i.e., 
^ms = 0- We separate Eq. (I11.3P into two kinds — the first kind with = that 
corresponds to exact non-trivial resonances, and the second kind that corresponds 
to non-trivial near-resonances and non-resonances. Since, in the summation, the 
first kind contains far fewer terms than the second kind, and all the terms are of the 
same order of magnitude, we will neglect the first kind in our analysis. Therefore, 
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Eq. flll.3p becomes 



N-l 



(11.5) 



l,m,s=l 



where the prime denotes the summation that neglects the exact non-trivial reso- 
nances. 

The problem of broadening of spectral peaks now becomes the study of the 
frequency spectrum of the dynamical variables bk (t) in thermal equilibrium. This is 
equivalent to study the two-point correlation in time of 



where the angular brackets denote the thermal average. Here we have used the 
Wiener-Khinchin theorem 



11.2 Wiener-Khinchin theorem 

Suppose we have a function of time f{t). Define its spectrum via the inverse 
Fourier transform 









Then, the Fourier transform takes the form 





Let us define the auto-correlation function 



Cit) = J f{Tyf{t + T)dT. 



(11.10) 



93 



Then the Wiener-Khinchin theorem states that 



(11.11) 



Proof: 



C(t) = y /(r)7(t + r)dr = ^ j f{u^)e-^-^^'^^^ r{u,)e^-^^du,du;,dT 
= 777-^ I fMe-'''''rM{2Ti)5{uj^-U2)dw^du2 



Although the proof of the theorem is quite simple, we can not use it for physically 
relevant signal, i.e., when f{t) is a stationary random process and has no Fourier 
transform. In this situation, we should define the autocorrelation function using 
expected value (over the Gibbs measure in our situation). However, for simplicity 
we will use the version that we have proven. 

11.3 Effective dynamics of the auto-correlation function C{t) 
Now we return to studying the dynamics of the correlation function of the 
interaction variable h{t). Under the dynamics (111.51) . time derivative of the two- 
point correlation becomes 




(11.12) 



QED. 



{hmm 




(11.13) 



where 
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In order to obtain a closed equation for Ck{t), we need to study the evolution of the 
fourth order correlator J^sit)- This correlator is a time-depended generalization 
of the similar correlator that was used in Chapter H] [Eq. (14 .Tp ]. We utilize the 
weak effective nonlinearity in Eq. (111.11) [22] as the small parameter in the following 
perturbation analysis and obtain a closure for Ck{t), similar to the traditional way of 
deriving kinetic equation, as in [5, 62] and which was discussed in Chapter |H We note 
that the effective interactions of renormalized waves can be weak, as we have shown 
in [22], even if the /5-FPU chain is in a strongly nonlinear regime. Our perturbation 
analysis is a multiple time-scale, statistical averaging method. Under the near- 
Gaussian assumption [Chapter H] , which is applicable for the weakly nonlinear wave 
fields in thermal equilibrium, for the four-point correlator, we obtain 

J^',(t)A^, = C,{t)Qm6y, + 5^51). (11.14) 

Combining Eqs. (111.21) and (111.141) . we find that the right-hand side of Eq. (111.131) 
vanishes because 

T'jisJtmms = ^- (11-15) 

Therefore, we need to proceed to the higher order contribution of J^si^)- Taking its 
time derivative yields 

j«,(t)Al = ([6r(t)6^(t)6.(t) + 6r(t)6^(t)6.(t) 

+ h:{t)h^{t%{t)]hm)^i,. (11.16) 

Considering the right-hand side of Eq. (111.161) term by term, for the first term, we 
have 



but)bsit)bm Ai- (11-17) 
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We can use the near-Gaussian assumption to split the correlator of the sixth order 
in Eq. flll.l7p into the product of three correlators of the second order, namely, 

{bi{o)bUtMt)b^{t)b;{t)b;)A'^^ = c,it)n^nj'^{6p; + 

Here, we have used that rim = Cm(0). Then, Eq. flll.l7p becomes 

{bKt)bUt)bsmm^'l = 4iTi^,C,(t)n^n,e-'^"=Aj^,. (11.18) 

Similarly, for the remaining two terms in Eq. flll.l6p . we have 

{b*{t)bUt)bsmmKs = -4if,rC,(t)n,n,e'^-^Al, (11.19) 

and 

{b*{t)bUt)bsmm^'l = -4iT;rC,(t)n,n„e^^«^A^,, (11.20) 

respectively. Combining Eqs. (I11.18p . (I11.19p . and (lll.2Up with Eq. (I11.16p . we 
obtain 

j^',(t)A^, = 4iT^',Cfc(t)e-'^"=*A^,(n^n.-n,n^-n,n.). (11.21) 

Equation flll.2ip can be solved for J^git) under the assumption that the term e~"^™»* 
oscillates much faster than Ck{t) [27]. We numerically verify [Fig. 111.41 below] the 
validity of this assumption of time-scale separation. Under this approximation, the 
solution of Eq. (111.210 becomes 

J^.^m'^s = 4T^',C,(t)Aj;^, ^^^rimns-ninm-nm,). (11.22) 

—UJ 

ms 

Plugging Eq. flll.22p into Eq. flll.lSp . we obtain the following equation for Ckit) 

at) = 8iCk{t)J2 [Tl) ^Is ki {nmn,-nins-ninm). (11.23) 
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Since in the thermal equihbrium rik is known, i.e., = (|&A;(t)P) = d/u}k [Eq- (17.831) ]. 
Eq. (111.231) becomes a closed equation for Ck{t). The solution of Eq. (111.231) yields 
the autocorrelation function Cfc(t) 



1 Ck{t) NAfcz 



;il.24) 



Using this observation, together with Eq. ( I11.2p . finally, we obtain for the thermal- 
ized /3-FPU chain 

Lm,s 



'msJ 



Equation (111.25^ gives a direct way of computing the correlation function of the 
renormalized waves a^, which, in turn, allows us to predict the spatiotemporal spec- 



trum |aA:(co')P. 



11.4 Analytical prediction vs numerical observation 

In Fig. 111.1( a). we plot the analytical prediction (via Eq. (111.251) ) of the spa- 
tiotemporal spectrum 

\dk{uj)\^ = Muj - Cok)\^ = d-'[Cit)]iuj - Uk). (11.26) 

By comparing this plot with the one presented in Fig. 111.1( b). in which the corre- 
sponding numerically measured spatiotemporal spectrum is shown, it can be seen 
that the analytical prediction of the frequency spectrum via Eq. (I11.25P is in good 
qualitative agreement with the numerically measured one. However, to obtain a 
more detailed comparison of the analytical prediction with the numerical observa- 
tion, we show, in Fig. I11.2[ the numerical frequency spectra of selected wave modes 
with the corresponding analytical predictions. It can be clearly observed that the 
agreement is rather good. One of the important characteristics of the frequency 
spectrum is the width of the spectrum. We compute the width W{f) of the spectrum 
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Figure 11.1: (a) Plot of the analytical prediction for the spatiotemporal spectrum 
|dA:(a;)p via Eq. (111.251) . (b) Plot of the numerically measured spa- 
tiotemporal spectrum |afc(ti;)p. The parameters in both plots were 
N = 256, /3 = 0.125, E = 100 and t] = 1.06, 6 = 0.401. r] and 6 were 
computed analytically via Gibbs measure. The darker gray scale cor- 
respond to larger values of |afc(a;)p in u-k space. [max{— 8, In |afc(u;)p} 
is plotted for clear presentation]. 



fiuj) by 



W{f) 



max^/(a;)' 



^11.27) 



In Fig. I11.3[ we compare the width, as a function of the wave number k, of the fre- 
quency peaks from the numerical observation with that obtained from the analytical 
predictions. We observe that, for weak nonlinearity (/5 = 0.125), the analytical pre- 
diction and the numerical observation are in excellent agreement. In the weakly 
nonlinear regime, this agreement can be attributed to the validity of (i) the near- 
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Figure 11.2: Temporal frequency spectrum |afe(co')p for k = 30 (left peak) and k = 50 
(right peak). The numerical spectrum is shown with pluses and the 
analytical prediction [via Eq. flll.25p ] is shown with solid line. The 
parameters were N = 256, (3 = 0.125, E = 100. 



Gaussian assumption, and (ii) the separation between the linear dispersion time 
scale and the time scale of the correlation Ck{t). This separation was used in deriv- 
ing the analytical prediction [Eq. (111.251) ]. However, when the nonlinearity becomes 
larger {(3 = 0.25 and j3 = 0.5), the discrepancy between the numerical measurements 
and the analytical prediction increases, as can be seen in Fig. 111.31 Nevertheless, 
it is important to emphasize that, even for very strong nonlinearity, our prediction 
is still qualitatively valid, as seen in Fig. (Ill.Sp . In order to find out the effect of 
the Umklapp scattering due to the finite size of the chain, we also computed the 
correlation [Eq. (111.251) ] with the "conventional" 5-function 6'^^ (i.e., without taking 
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Figure 11.3: Frequency peak width Ty(|afc(ti;)p) as a function of the wave number 
k. The analytical prediction via Eq. (I11.25P is shown with a dashed 
line and the numerical observation is plotted with solid circles. The 
parameters were = 256, E = 100. The upper thick lines correspond 
to jS = 0.5, the middle fine lines correspond to j3 = 0.25, and the lower 
solid circle and dashed line (almost overlap) correspond to /? = 0.125. 



into account the Umklapp processes) instead of our "periodic" delta function A^^. 
It turns out that the correlation time is approximately 30% larger if it is computed 
without Umklapp processes taken into account for the case = 256, P = 0.5, 
E = 100. It demonstrates that the influence of the non-trivial Umklapp resonances 
is important and should be considered when one describes the dynamics of the finite 
length chain of particles. Finally, in Fig. 111.41 we verify the time scale separation 
assumption used in our derivation, i.e., the correlation time of the wave mode k is 
sufficiently larger than the corresponding linear dispersion period = 27T/ujk. In the 
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case of small nonlinearity {(3 = 0.125), the two-point correlation changes over much 
slower time scale than the corresponding linear oscillations — the correlation time 
is nearly two orders of magnitude larger than the corresponding linear oscillations 
for weak nonlinearity /3 — 0.125, and nearly one order of magnitude larger than the 
corresponding linear oscillations for stronger nonlinearity P — 0.25 and P — 0.5. 
This demonstrates that the renormalized waves have long lifetimes, i.e., they are 
coherent over time-scales that are much longer than their oscillation time-scales. 
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Figure 11.4: Ratio, as a function of k, of the correlation time of the mode k 
to the corresponding hnear period tk = 271 /uk- Circles, squares, and 
diamonds represent the analytical prediction for /5 = 0.5, P = 0.25, 
and P = 0.125 respectively. Solid circles, pentagrams, and triangles 
correspond to the numerical observation for P = 0.5, (3 = 0.25, and 
(3 = 0.125 respectively. The parameters were N = 256, E = 100. The 
ratio is sufficiently large for all wave numbers k even for relatively large 
f3 = 0.5, which validates the time-scale separation assumption used in 
deriving Eq. (111.221) . The comparison also suggests that for smaller (3 
the analytical prediction should be closer to the numerical observation, 
as is confirmed in Fig. 111.31 



CHAPTER 12 
Discrete Breathers in /3-FPU chains 



We have studied the interactions of renormahzed waves in the FPU systems in 
thermal equihbrium. As we have seen, under any strength of nonhnearity, the 
renormahzed waves in thermal equilibrium can be regarded as quasiparticles in the 
wave-number space. The question arises: are there any localized excitations of the 
FPU lattice in the physical space. In the last decade, discrete breathers (DB) as spa- 
tially localized, time periodic lattice excitations were discovered [19]. Arising from 
energy localization in nonlinear lattices, they play important roles in many dynamics 
in fiber optics, condensed matter physics and molecular biology [17]. The existence 
of DBs has been addressed rigorously [4]. Important conceptual issues naturally 
appear, such as what is the role of DBs on the route to equilibrium [16] and how 
do they manifest in thermalization of the FPU system? Resolution of these issues 
will certainly provide deep insight into the fundamental understanding of route to 
thermalization for general nonlinear physical systems. Most of the results regarding 
DBs in (3-FPU chains have so far only addressed their behavior in the transient state 
of weakly nonlinear regimes before thermalization occurs [15,54]. Here we present 
numerical evidence of DBs in the FPU chains in thermal equilibrium. However, 
before we do that we will give a short introduction to the large subject of breather 
solutions in both continuous and discrete models. A more complete review of the 
topics can be found in [19]. 

12.1 General introduction to breathers 

Breathers are spatially localized, time-periodic solutions of a PDE or of a 
lattice system. Spatial localization usually means exponential decay of the ampli- 
tude of a breather from its center. A well known example of a PDE that possesses 
breather solutions is a particular case of Klein-Gordon equation 

At = Ci:,,-F{i:), (12.1) 
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with the choice F{iIj) = sm{ilj) and which is referred to as sine-Gordon equation. 
The breather solution for the sine-Gordon equation is given by 

^, = 4tan-i— (12.2) 
ujch[mx) 



where u = \f\^^m^ . However, these breather solutions are structurally unstable 
in the following sense. For most of the perturbations of the nonlinear term -F(V'), 
the breather solutions do not survive. For instance, it was proved by Segur and 
Kruskal [55] that for the choice = — -|- (a so called 0^ system) breather 

solutions do not exist. The intuitive underlying reason for breather solutions of 
PDEs to be structurally unstable arises from the fact that the linear spectrum of 
continuous PDEs is unbounded and the frequencies of the breathers would most 
probably be in resonance with the linear spectrum. In particular, the linear disper- 
sion relation of the Klein-Gordon model is [19] 



^/Ck^ + F'{ij = 0). (12.3) 



This dispersion relation is unbounded and, therefore, it is difficult to find breather 
solutions in general PDEs except for a comparably few non-generic cases. 

The situation changes drastically when we turn to the discrete nonlinear mod- 
els. Consider a discretization of the Klein-Gordon model 

iij + 6^(2U„ - Un+l - Un-l) - f{Un) = 0. (12.4) 

Here the linear dispersion relation reads [19] 

ujk = Vl + 2&2(l-cos(A;)). (12.5) 

We note, that unlike the spectrum (112. 3p of the continuous model, the spectrum (112. 5p 
of the discrete model is bounded from both below and above, which allows for the 
existence of breathers if their frequencies lie outside the linear band. It was demon- 
strated in [41] that systems like (112. 4p have discrete breather solutions. The result 
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Figure 12.1: Spatial localization and periodicity in time (usually with oscilla- 
tions of the alternating wings) are characteristic properties of discrete 
breathers. 

proved in [41] is the most general known so far. Indeed, it was actually proved that 
almost any nonlinear lattice possesses discrete breather solutions if it satisfies some 
non-resonance conditions. 

In Fig. 112.11 we can see an idealized form of the discrete breather solution 
of a nonlinear lattice. The characteristic properties of the DBs are: (i) spatial 
localization — in practice, highly localized DBs usually involve just 4-5 sites of the 
lattice — and (ii) time-periodicity — the subsequent sites oscillate in an alternating 
manner with the frequency, higher than the linear band of the spectrum. The 
situation when the frequency of the breather lies below the linear band is also 
possible, however, we will not investigate this kind of breathers in this thesis since 
the linear dispersion of the FPU chains satisfies < ujk < 2. 
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12.2 Discrete Breathers in the /5-FPU chain in transient to 
the thermal equihbrium 

The formation of Discrete Breathers in transient to thermal equihbrium in the 
weakly nonlinear /3-FPU chain is studied in detail in [15]. Here, we provide some 
of the results of this work and in the next Section we will discuss the existence of 
DB's in the /5-FPU system after it has reached the thermal equilibrium state. 

As opposed to the original FPU experiment, the authors of [15] considered the 
initial condition with the highest Fourier mode excited. It is shown that this mode 
is modulationally unstable and gives rise to the energy localization in the form of 
DBs. In order to observe these DBs numerically, one can integrate the dynamical 
equations of motion (18. 2p with the initial condition of a zig-zag form 



where a is the amplitude of the initial condition. This initial condition is usually 
referred to as vr-mode. vr-mode is an exact solution of the /5-FPU lattice [49] and, 
therefore, in order to destabilize it, small noise needs to be introduced to the initial 
condition — the random perturbation of the order 10~^^ was used in Pj{0)- 

It is convenient to detect energy localization using the function L{t) given in 
Eq. (18. 9p . Recall that L{t) is of the order 1 if energy is nearly uniformly distributed 
along the chain and of the order if energy is concentrated around a few sites of 
the chain. In Fig. 112.21 we show the time evolution of the chain qj{t) (panel (a)), 
together with the time evolution of L{t) (panel (b)). First the oscillations of the 
chain follow the initially given zig-zag pattern, which corresponds to the uniform 
energy distribution with L{t) = 1. Then the random initial perturbation induces 
modulation instability and energy starts to localize. At that time we can observe 
the development of DB's (dark strips in Fig. 112.2( a)) and the increase of the value 
of L(t) (Fig. 112. 2f b)). Then the localized states disperse their energy and gradually 
disappear, and the system approaches the energy equipartition state. At that point 
L{t) goes to some steady state value. In Fig. 112. 3[ we demonstrate the snapshots of 



0, 



i-iya, 



(12.6) 
(12.7) 
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the chain, i.e., q{j), at the three different stages of the chain evolution that we have 
just described. In Fig. 112. 3T a). the initial condition is shown. In Fig. 112.3T b). we see 
the energy localization state — DBs are developed (enclosed in the red rectangle). 
Note the amplitude of qj. And finally, in Fig. 112.3( c). the snapshot of the chain at 
the thermalized state is shown, when the chain only consists of the renormalized 
linear waves. 

It is known [19] that spatially localized lattice solutions have frequencies that 
are outside the linear band of the system. Otherwise, they would resonate with the 
spatially extended linear waves. The linear band of the /5-FPU chain is given by the 
dispersion relation Uk = 2sin{TTk/N). Therefore, frequencies of the DBs must be 
greater than 2 as we have discussed in Section [12.1[ We confirm this fact numerically. 
In Fig. 112.41 we plot the spatiotemporal spectrum of the /?-FPU chain at the time 
period when the DBs are developed (corresponds to t G [2000,8000] in Fig. 112. 2p . 
It can be seen in Fig. 112.41 that the spatiotemporal spectrum essentially consists 
of two parts: (i) the linear band given by Uk (we also notice renormalization due 
to the nonlinearity, see Chapter [7] for details) and (ii) the higher frequency region 
(enclosed in the black rectangle) that corresponds to DBs. In order to compute the 
spatiotemporal spectrum, we have used the time series with a length of 16384 time 
steps, and each time step was 0.1 time units. 

Now that we have observed the existence of DBs in the transient to thermal 
equilibrium, we turn to investigating the existence of DBs in the /3-FPU system 
after it has reached the thermal equilibrium state. 




Figure 12.2: qj{t) as a function of j and t (panel (a)) and energy localization L(t) 
(panel (b)) computed via Eq. fl8.9p . The chain was modeled for N = 
128, (3 = 0.1, and a = 0.8 {E/N ~ 1.44). 
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(a) Modulational instability t=10 
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(b) Energy localization t=2000 




20 40 60 80 100 120 

(c) Energy equipartition t=10^ 
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Figure 12.3: Panel (a): initial state. Panel (b): energy localization state. The 

breather is enclosed in a red rectangle. Panel (c): energy equipartition 

state. The chain was modeled for = 128, /5 = 0.1, and a = 0.8 
(E/N ~ 1.44). 



109 




UJ 



Figure 12.4: Spatiotemporal spectrum demonstrates the oscillations with frequen- 
cies outside the linear band (region enclosed in a black rectangle is 
separated from the sin-like linear band.) The chain was modeled for 
= 128, /3 = 0.1, and a = 0.8 {E/N ~ 1.44) and the spectrum was 
computed via the time series of the length 16384 time steps with each 
time step of the length 0.1 time units. 
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12.3 Discrete Breathers in the thermal equihbrium 

As we have already seen in Chapter [71 the thermahzed state of the /5-FPU 
chain is characterized by the existence of the renormahzed waves. Moreover, in 
Section 112. 2[ we have discussed the existence of DB excitations that were observed 
during transient stages towards thermahzation [15]. Here we show via numerical sim- 
ulation that DBs actually persist and coexist with renormahzed waves in the ther- 
mahzed state. Thus, in the thermahzed /3-FPU, there are two kinds of quasi-particle 
excitations, one localized in fc-space as renormahzed nonlinear waves/phonons, and 
the other, localized in x-space as DBs. 

As nonlinearity increases, the duration of transient becomes shorter. After 
the energy redistributes among all the modes to achieve thermal equilibration, our 
simulations show that the spatially localized, high frequency excitations still exist. 
These DBs can interact with each other and may be destroyed by collision processes 
with other DBs or with the renormahzed waves. The spatial structure of these 
excitations very much resembles the idealized breather oscillations (Fig. 112.11) in the 
absence of spatially extended waves: they "live" above the high frequency edge of 
the dispersion band and their lifetime is sufficiently long (on the order of 10-100 DB 
oscillations). Therefore, they behave like a quasiparticle. Note that, under certain 
conditions, supersonic solitons may arise from the /5-FPU system as another kind of 
localized excitations [32,64]. However, they were not observed in our thermahzed 
system. 

Figure 112.51 is the energy density plot which shows the time evolution of en- 
ergy of each particle for the transient (recording starting time Tq = 5 x 10^) and 
thermahzed (Tq = 5 x 10^) states, respectively. Fig. 112.5( c) and (d) display the 
energy as a function of site at Tq corresponding to Fig. 112.5( a) and (b) respectively. 
In the transient case (Fig. 112.5( a)) the spatially localized objects (dark stripes) that 
carry sufficiently large amount of energy are clearly observed. Fig. 112.5( c) is a snap- 
shot of the energy density plot (112. 5( a)) at Tq. Here the DBs are seen as localized 
peaks [15]. After thermahzation the spatial structure looks different (Fig. 112.5( b)). 
The system now consists of the renormahzed waves (straight cross-hatch traces in 
Fig. 112.5( b)). On the top of these waves, the localized structures similar to DBs 
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manifest themselves as the wavy dark trajectories (in Fig. 112.5( b)). Although the 
snapshot (Fig. 112.5( d)) of the energy density plot (Fig. 112.5( b)) indicates that in 
thermal equilibrium the energy is more evenly distributed among particles, spatially 
localized structures are clearly observed. 

Since there are renormalized waves in the system, which also carry energy, 
we need to find a way to distinguish between these waves and DBs. We use a 
frequency filter that cuts out the lower side of the Fourier spectrum and leaves 
the high frequency part unmodified, i.e., f{g{n,t)) =ReF~^{H^{F{g))), where Re 
denotes the real part, F is a time Fourier transform, eliminates all frequencies 
below uJcut and g{n,t) is a dynamical variable that is being filtered. By applying 
this filter to the displacement g„ to obtain ql = f{qn), we can show the existence 
of DBs even for strong nonlinearities, for example, (3 = 25. Fig. 112.6( a) shows a 
clear example of a DB excitation reconstructed using the filtered with Ucut = 
7. Fig. 112.6( b) shows a typical spatial profile of the DB taken from Fig. 112.6( a). 
which strongly resembles the idealized DB [7]. Finally, in Fig. 112. 7[ we present the 
evidence that there is a turbulence of DBs, which chaotically ride on renormalized 
waves. The corresponding energy density distribution along with the distribution of 
filtered displacement in a zoomed region qlit) is displayed in Fig. 112. 7T a) and (b), 
respectively. After the lower modes from the displacement g„ are filtered, one can 
clearly observe that the remaining high frequency oscillations are spatially highly 
localized, with the same characteristics as an idealized breather. The detailed time 
dynamics of the DB shows the main characteristics of breathers: the values of 
change signs periodically (as indicated by the alternating white and black spots 
along the trajectory) as the DB moves in space, with a spatial span of 2 or 3 sites 
only, as seen in Fig. 112.7( b). 

To summarize, we have demonstrated numerically, that DBs persist in the 
thermal equilibrium state of the /?-FPU chains. However, in order to detect them, 
the linear dynamics of the chain has to be filtered out. The remaining excitations 
appear to be localized in space and periodic in time, i.e., they satisfy the definition 
of DB. 



112 





Figure 12.5: Energy density evolution of the transient state (a) and (c), of the ther- 
mahzed state (b) and (d), respectively (/3 = 1 and H = 200). In (a) 
and (b) the darker strips correspond to high energy localizations, (c) 
and (d) are the snapshots of energy density. 
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Figure 12.6: (a) Evolution of a discrete breather in thermal equilibrium, (b) typical 
snapshot of the breather. f3 = 25, H = 200. 




Figure 12.7: Turbulence of discrete breathers (N=1024): (a) evolution of energy 
density, (b) zoomed in ql{t) of the area indicated by the rectangle in 



CHAPTER 13 
Conclusions 



In this thesis, we have investigated one of the most famous problems of the nonhn- 
ear science, the celebrated FPU system. We have studied the problem from both 
dynamical and statistical angles. 

In the beginning of the thesis, we have given some historic notes about the FPU 
problem as well as an overview of the topics that are widely used in treating nonlinear 
problems: classical mechanics, statistical mechanics, wave turbulence, numerical 
methods, and chaos. 

Then, we turned to the FPU problem and studied the statistical behavior in 
thermal equilibrium. We have extended the notion of normal modes to the nonlinear 
system by showing that regardless of the strength of nonlinearity, the system in ther- 
mal equilibrium can still be effectively characterized by a complete set of renormal- 
ized waves, in the sense that those renormalized waves possess the Rayleigh- Jeans 
distribution and vanishing correlations between different wave modes. In addition, 
we have studied the property of dispersion relation of the renormalized waves. The 
results we obtained in Chapter [7| are general and can be applied to the large class 
of nonlinear systems with the nearest-neighbor interactions in thermal equilibrium. 

We have further focused our attention on the /3-FPU chain, which is character- 
ized by the fourth order potential. We have confirmed that the general renormaliza- 
tion framework that we discussed above is consistent with numerical observations. 
In particular, we have shown that the renormalized dispersion of the thermalized 
/3-FPU chain is in excellent agreement with the numerical one for a wide range 
of the nonlinearity strength. We have further demonstrated that the renormalized 
dispersion is a direct consequence of the trivial resonant interactions of the renor- 
malized waves. Using a self-consistency argument, we have found an approximation 
of the renormalization factor via a mean-field approximation. In addition, we have 
used the multiple time-scale, statistical averaging method to obtain the theoretical 
prediction of the spatiotemporal spectrum and demonstrated that the renormalized 
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waves have long lifetimes. 

Moreover, we studied the FPU system from the particle interaction point of 
view. In particular, we have investigated the existence of the discrete breather 
excitations in the /9-FPU chain. We have numerically demonstrated that the discrete 
breather solutions that were observed previously in the transient to the thermal 
equilibrium still persist even when the system reaches thermal equilibrium. 
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APPENDIX A 
KdV equation as a continuous approximation of the FPU 

chains 

In this Appendix, we derive the KdV equation as a continuous approximation of the 
a-FPU chains [30]. Consider an a-FPU chain, which is given by a Hamiltonian 

N 2 

where the potential is of the form 

Pis) = ^s^ + %^ (A.2) 
Then the equation of motion becomes 

mqj = K{qj+i - 2qj + qj_i) + a ((g^+i - q^f - (g^ - qj^if) ■ (A.3) 
Let us rescale time t and the nonhnearity parameter a via 

t ^ J^t, (A.4) 
V m 

a — > —a. (A. 5) 

n 

Now, Eq. (lA.Sp becomes 

(ii = - "^Ij + + " - (Ijf - i<lj - (Ij-if) ■ (A-6) 

Next, we denote Uj = qj — qj-i and rewrite Eq. (]A.6P as 



y, = - 2F{y,) + (A.7) 
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where 



F(s) = s + as' 



(A. 



Suppose the chain has length L and spacing h. We consider the continuous hmit 
N (yo and /i — > such that Nh = const = L. Denote x = nh and y{x) = Un- 
Then, y{x) can be regarded as a function of the continuous-valued variable x. Using 
the Taylor expansion, we have 

Ffe.O = F(yinh + h)) = FU.)) + h?I^ + \h-?^I^ + ... 



Using this formal notation, Eq. ( ]A.7I) can be rewritten as 



(A.9) 



(A.IO) 



or, equivalent ly 



=4sinh^ ( % ] F{y{x)) 



(A.ll) 



Next, we use the Taylor expansion of the sinh function up to 0{h^) and neglect 
0{ah'^) and 0{h^). Here we use a ^ 1. Equation flA.lip becomes 



y{x) = 4 



3 V 2 



{y{x) + ay{xf). 



(A.12) 



After applying the differential operators, Eq. flA.12p takes the form 



y = h y^^ + ah {y^)xx + —yxxxx- 



(A.13) 
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Let us rescale time, space and displacement variables via 



t Vm, (A.14) 



Vl2 

X — > — — X, (A. 15) 

n 

y ay. (A. 16) 



Then, Eq. (1A.13P becomes 

■ij = (y + + yxx)xx- (A.i7) 

Now, let is return back to the initial variables m„. Similarly to y{x), we denote u{x) 
to be a function of a continuous variable x 

u{x) = u{nh) = Un- (A. 18) 

Then, using the same rescaling u —>■ au, we obtain the following connection between 
y and u 

y = VT2u^. (A. 19) 

Combining Eqs. ( 1A.17I) and ( lA.19p . we obtain the Boussinesq equation 



ii = (1 + 2ux)uxx + Uxxxx- (A.20) 
Let us first consider the case without the dispersion. We can rewrite Eq. ( ]A.20[) as 

u - H\u,)u^^ = 0, (A.21) 



where H{s) = a/I + 2s. Next, we introduce new variables 



w = Ux, 

(A.22) 

V = Ut. 
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Then Eq. (IA.2ip becomes a system 

Wt-Vx =0 



(A.23) 
vt - H'^{w)w^ = 0. 



Let us make a linear transformation of Eq. flA.23P in the following way. First, we 



add the first equation multiplied by H with the second one. Secondly, we subtract 
the second equation from the first multiplied by H. Thus, we obtain 



Hwf — Hvx + Vt — H'^Wx = 0, 
—Vt + Hwt — Hvx + H'^Wx = 0. 



System (]A.24p is equivalent to the following system 




where 

s = V + / H(r])dr], 
Jo 

r = —V + / H(r])dr]. 







Therefore, we have 



(A.24) 



(A.25) 



Let us express w in terms of r and s 

r + s = 2 / H{7])dr] = 2G{w). (A.26) 
Jo 



w = G-'( ^ ) . (A.27) 
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Then, the second equation in (1A.25P becomes 



rt + H{G-'r-±^] I =0 (A.28) 



Now, we consider the case with the dispersion. The same procedure gives us 

■^t H Sx Wxxxi 

(A.29) 

Tt + HTx = —Wxxxi 

The LHS of the first equation in (1A.29P is a complete derivative along the charac- 
teristics 

ds 

= uJxxx, (A. 30) 

along 

^ = -H. (A.31) 

dt ^ ' 

in Eq. flA.301) . we assume that the RHS is small for long waves. Therefore, s is a 
constant. Without loss of generality, we assume that s = 0. Then, we have 

1 r -I- s 

^l^d^ = 3 ((1 + 2«^)'^' - 1) = 

Taking into account that s = 0, we find 



Therefore, we obtain 



3r ^ 



F{w) = ( Y + 1 ) ■ (A.34) 
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For small u we have the following approximations 



r 

^ = 2' 
r 

H = - + 



(A.35) 



Substituting Eqs. (1A.35P into the second equation in (IA.29p . we have 



n + + l) + ^r^xx = 0. (A.36) 

And finally, making the following changes of variables x ^ x — t and then t — >■ t/2, 
we obtain the KdV equation 

rt + rr^ + r^^^ = 0. (A.37) 

Therefore, we have shown that the KdV equation is indeed a continuous approxi- 
mation of the a-FPU chains in the small amplitude and long wavelength regime. 
Similarly, one can show that the so called modified KdV equation 

n + r^r, + r,,, = (A.38) 



is a continuous approximation of the /3-FPU chain. 



APPENDIX B 

Computation of the Lyapunov exponent of the /3-FPU chain 



In order to investigate the chaotic structure of the /?-FPU chain given by Eq. flS.ip . 
we measure the Lyapunov exponent. We consider small perturbations 6j and ej of 
the dynamical variables qi and pi, respectively. As we have studied in Section [6731 
a system exhibits chaotic behavior if the small perturbations of the dynamical vari- 
ables grow exponentially with time. In order to describe the dynamical behavior of 
the perturbations 6j and ej, we linearize Eqs. (18. 2 p and then numerically study their 
evolution 

~ ' 

Bj = {6j+i - 2Sj + Sj_i) - 3(3{{qj - qj+iYidj - 6j+i) + {qj - qj-if{Sj - Sj^i)). 

(B.l) 

Note that we have to solve Eqs. (18. 2p together with Eqs. flB.lP . The procedure of 
computing the Lyapunov exponent is the following. We take some initial condition 
{qi{0) , pi{0)) and initial perturbation (ej(0), 5j(0)). We choose the perturbation to 
be very small, e.g. its I2 norm (which we denote as do) is of the order 10'^°. Then, 
let the dynamical variables pj and qj and the perturbations 6j and Ej evolve for n 
time units according to Eqs. (18. 2p and Eqs. (IB.ip . respectively. After n time units, 
we measure the norm of the perturbation again (denote it as di). Let us call 

hi = -\og^. (B.2) 

n do 

Now, the perturbation has to be rescaled to make it have the norm do again. The 
direction of the perturbation vector {e, 6) has to stay the same, as it was before 
rescaling, only the absolute value takes its initial value. After we repeat this proce- 
dure m times we obtain the values of h^"^^ = ^ hs- The Lyapunov exponent 
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can be then estimated as 

h= lim = lim — y^hj. (B.3) 

Note that the forth order Runge-Kutta method can be used for solving Eqs. fIB.ip . 

Therefore, we have presented a numerical algorithm of computing the Lya- 
punov exponent of the /3-FPU chain. 



